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by 
Norbert H. Doerry 
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the requirements for the degree of Doctor of Philosophy in the field of Naval Electrical 
Power Systems. 
ABSTRACT 

Naval shipboard electric power systems are transitioning from the relatively simple 
distribution of ship service electric power to extremely complex integrated electric drive 
(IED) systems. The optimal design of warships employing IED is presently hampered by the 
lack of existing simulation computer tools for analyzing the highly coupled and controlled 
electro-mechanical systems characteristic of IED. As a first step in the development of a 
viable computer simulation tool, the numerical algorithm testing program WAVESIM was 
created. 

The key contnbutions of WAVESIM are the systematic treatment of waveforms as an 
abstract data type, the development of the terminal description of devices, and the use of 
structural jacobians in system reduction. 

WAVESIM represents variables by waveforms consisting of a vector of coefficients 
and a waveform type code indicating how the coefficients should be interpreted. The 
principal advantage of using waveforms over conventional discrete point methods is the 
avoidance of unstable integration techniques since for most waveform types, integration and 
differentiation are linear matrix operations. 

Devices are described in WAVESIM by relationships between terminal interface 
variables. WAWVESIM recognizes two types of terminals: normal terminals having both 
potential and flow variables, and information terminals having only a potential variable. In 
this manner, WAVESIM can simulate processes involving both energy transfer and control 
Signals. 

WAVESIM extends the structural jacobian matrix concept to reflect the properties of 
the dependence of system equations on system variables. The system structural jacobian 
matrix is constructed from the constitutive device structural jacobian matrices and is used to 
identfiy a sequence of smaller blocks when can be solved consecutively for all the system 
variables. 

To demonstrate and verify the capabilities of WAVESIM, several simulations were 
conducted. In all simulations, WAVESIM provide results matching data generated by other 
simulation methods. 
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Chapter 1 Introduction 


A revolution is occurring in modern warship design. The conventional mechanical 
transmission train for transferring power from the prime movers to the ships screws will be 
replaced in future warship designs by an integrated electric drive (IED) system. While 
electric drive is not a new concept, the IED approach differs significantly from previous 
electric drive implementations in that both propulsion power and ship service power (60 HZ 
440 Volts AC) are derived from the same prime movers. The resulting flexibility in the 
arrangability of the ship, in the design of the electric plant, and in the power available to 
combat systems provides the naval architect with many opportunities for significantly 


improving the combat effectiveness of modern warships. 


Designing a ship taking full advantage of the opportunities afforded by IED is not an 
easy task or even obvious. The ship designer has no precedent for an IED ship let alone the 
design of an IED electric plant. Instead, the designer must rely heavily on simulations of 
proposed systems to evaluate the soundness of the design. For the electrical power system, a 


suitable simulation environment must be capable of addressing these questions: 


Will the Electric Power System Work? 

This is the ultumate question which needs to be answered. Unfortunately defining 
the term work is not an easy task, nor is assuring a system will work under all operating 
conditions. A strict time domain simulation only provides a solution for a given set of 
operating conditions. Generalizing the results of relatively few simulations to all 
operating conditions is both necessary and prone to catastrophic failure. Hence more 


than just a time response is usually needed. 


How Will the System React to Major Disturbances and Faults? 

The primary design goal for shipboard electric power systems is continuity of 
power. To this end, the response of the system to abnormal events such as grounds, 
stalled motors, and inadvertent opening of breakers is crucial to evaluating the success 


of the electric power system design. 


How Will the System React to Severe Dynamic Conditions? 
A number of normal events can cause severe dynamic responses within the 
system. Rapidly changing the propulsion motor speed or direction, discharging pulse 


weapons, or Starting large motors are all examples of normal dynamic events. 
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Is the System Stable During a Given Dynamic Scenario? 
One import aspect of a system that works is its stability. The system should 
remain stable during all normal dynamic conditions and for as many disturbances and 


faults as possible. 


What is the Stability Margin? 
Some measure of how stable the system is should be provided to assist in 


generalizing the findings of stability about one scenario to other related scenarios. 


What is the Sensitivity of the Results to Parameters? 

The generation of models for a dynamic system simulation requires some 
estimation of device parameters. Knowledge of the sensitivity of the simulation results 
to parameter estimation errors is crucial for correlating simulation results with what can 


be expected from the physical system. 


How Correct are the Answers Provided to the Above Questions? 

Simulations generally use numerical methods to generate solutions. Careful 
control of error propagation is very unportant in ensuring accurate conclusions can be 
drawn from the simulation results. Some form of feedback should be provided to the 


operator as to the confidence level of the results. 


These requirements for performing time domain simulations of proposed and existing 
electric power systems found on United States naval warships can be quite challenging. The 
size, complexity, and strong coupling of components all conspire to make the simulator’s 
task difficult. At first glance, one would think the simulation programs designed for the 
commercial power utilities would be sufficient for handling the smaller shipboard systems. 
Unfortunately, this is not the case due to the following differences of the shipboard system 


from commercial power systems: 


Variable Frequency 

Frequency cannot be assumed constant. Many IED designs have the generators, 
motors, and ship service power all operating at different frequencies to optimize the 
performance of individual components. Frequency changers are employed to convert 
the power from one frequency to another. Even the ship service system onboard 
mechanical drive ships can experience frequency fluctuations lasting up to 2 seconds. 
The limited rotational intertia of the prime movers and generators allows for rapid 


accelerations and decelerations of the shaft and corresponding frequency fluctuations. 
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Lack of Time Scale Separation 

The principal time constants of controls, machine dynamics, and electric 
dynamics all fall within the same general range of milliseconds to seconds. The 
practice of decomposing the problem by time scale separation often used in analyzing 


commercial power systems becomes much more difficult. 


Load Sharing instead of Power Scheduling 

The commercial power utilities operate by scheduling the power delivered by 
each of the generating units. The mismatch between scheduled power generation and 
the actual load is met by a swing generator. Onboard ship however, both real and 
reactive power are shared equally among all paralleled generators through the very fast 
exchange of load sharing information. This fast exchange of information strongly 


couples the dynamics of all the paralleled generators. 


Short Electrical Distances 

The distances onboard ship are so short (under 1000 ft) as to make the modelling 
of transmission lines unnecessary for most simulations and to trivialize the load flow 
problem which is so important to the commercial power sector. The short electrical 
distances also strengthen the coupling of the various subsystems making up the 


electrical power system. 


Load Dynamics 

Commercial utilities usually assume loads are either consuming constant real and 
reactive power, or are Constant impedances. Shipboard systems however, must account 
for dynamics of loads such as propulsion motors, large pumps, pulsed loads, propeller 
dynamics, and ship dynamics. Furthermore, the supervisory level control envisioned 
for future designs may have the ability to control aspects of the loads in addition to 


generation. 


Tighter Control 

Because a ship is relatively small, a higher level of control can be exercised over 
the shipboard power system than can be exercised in the commercial power industry. 
This higher level of control strengthens the dynamic coupling of components of the 


system and complicates simulation efforts. 


Clearly, shipboard systems are considerably different from commercial power systems; 
and the inapplicability to shipboard power systems of simulation techniques developed for 


commercial systems should not be surprising. Other simulation tools exist but for one or 
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more reasons, all are ill-suited for simulating shipboard systems. A review of currently 
available commercial software for solving lumped parameter systems reveals no program 
currenty exists suitable for fulfilling the needs of simulating shipboard electric power 


Systems. 


Circuit Simulators 

As will be discussed in following sections, circuit simulators almost universally 
describe devices in terms of branch voltages and curents. The constitutive 
relationships are based only on the relative difference of the terminal variables and can 
not depend on the actual nodal potentials. Furthermore, the flow variables must be 
conserved on the device level. While these restrictions are not of concer for electrical 
networks, they are a bit constraining on electro-mechanical systems where one would 
like to deal with energy transformations in a device by employing equations which do 
not conserve the flow variable on the device level. The torque on the input shaft of a 
gearbox for example, does not equal the torque on the output shaft. Even electric 
power models where the flow variable is power and the potential variable is voltage can 
best be described by constitutive equation which do not enforce conserving power by 
ignoring the power converted to heat through resistive losses. 

Many circuit simulator also have problems modelling the transfer of information 
which is common in systems employing control systems. Information has only 


potentials and no flows associated with it. 


Signal Analysis Software 

Signal Analysis Software is often used to simulate control systems but often have 
difficulty simulating energy transfer. In particular, these programs often are incapable 
of solving implicit equations which are typically created by writing Kirchhoff’s Current 
Law when developing systems. Instead much effort must be expended to ensure the 


models have the appropriate input and output variables for a given system to be built. 


Commercial Power Utility Analysis Programs 

Software for analyzing commercial power universally do not apply to shipboard 
Systems due to several key differences. The lack of transmission lines, rotational 
inertia, time scale separation of dynamics and the presence of tightly coupled control 
loops are all features of the shipboard system which prevent the use of the commercial 


power system analysis techniques [5] [9] [10] [11]. 
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General Differential Equation solvers 

Most general differential equation solving programs cannot handle implicit 
equations very well. While the development and interconnecting of models into 
systems is possible, the process can be quite cumbersome [12] [13]. Dynamically stiff 


systems also pose serious challenges to the general differential equation solvers. 


Hybrid Computers 

Hybrid computers for studying electrical power systems can answer many of the 
desired questions for small shipboard systems. Unfortunately, hybrid computers are 
very expensive and limited in the size of problems which can be addressed. Presently 
the only hybrid computer in the United States suitable for shipboard power system 
studies is located at Purdue University. While this machine 1s capable, the needs of the 
IED program will require digital computer techniques for performing the desired 
studies. [92] [93] [94] [95] [96] 


As part of an effort to fill the need for simulating shipboard power systems, the 
WAVESIM program was specially created to develop applicable simulation techniques. 
Before discussing the numerical methods employed in WAVESIM, a review of existing 


methods is appropmiate. 
1.1 Simulation Process 


The process of simulating a physical system can be broken into three steps. First, a 
system of equations descnbing the component device constitutive relationships as well as 
the network constraints must be developed. While the network constraints are purely linear 
algebraic equations, the constitutive equations can be nonlinear and dynamic. Together, a 
system of nonlinear differential algebraic equations is generated. The next step is the 
conversion of the system of differential algebraic equations into a sequence of purely 
algebraic equations. Common integration techniques include the forward and backward 
Euler methods, Trapezoidal rule integration, and the Runge-Kutta methods. Finally, the 
nonlinear algebraic system 1s solved either through the Newton-Raphson method or through 


one of several relaxation techniques. 


Before describing several methods for generating and solving the system of equations 
corresponding to a physical system, the difference between the branch description and 
terminal description of devices should be detailed. The branch description of devices 
requires all the constitutive relationships be based on the relative difference between 


terminal potentials and all flows entering a device also leave the device. Hence for a two 
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terminal device, there is a single branch potential variable and a single branch flow variable 
associated with it. In the terminal description, the potential and flow associated with each 
terminal are variables. A two terminal device would then have four variables associated 
with it: two flow variables and two potentials. The terminal description allows the 
constitutive equations be a function of the actual values of the terminal potentials and not 
only of the potential difference. In other words, the potential reference can be set at the 
system level and not necessarily on the device level. Furthermore, the flows are not 
required to be conserved. A gear box for example, has differing torques entering and 
exiting it. The branch description method requires a four terminal model of the gear box 
while the terminal description requires only two terminals. In either case the result would 
be four variables descnbing the gearbox, but the branch description requires an explicit 
declaration of the device potential reference while the terminal description uses an implicit 


system wide potential reference. 


Branch Description vs. Terminal Description 
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1.2 Developing System Equations 


The normal method of describing a dynamic system is to generate a consistent set of 
possibly nonlinear differential algebraic equations and arrange them into the following 


canonical form: 
Ci =f(%,¥,u) 
O=g(x,y,u) 
where x is the vector of dynamic state variables, y is the vector of state variables with no 


associated dynamics staes, and u is the vector of system inputs. This system of differential 
algebraic equations (DAE) can be generated several different ways with the most common 
being the Sparse Tableu, Modified Nodal Analysis, and the standard load flow method. 


The method employed in WAVESIM does not extract the differential equations from 
the device constitutive equations but instead forms a system of algebraic equations of the 


form: 
O= bie Ge aiean) 


where x is the vector of the system variables and g,() is a device function having subsets x; 
and 
z, 


of x and u as arguments. The functions g,{) have the dynamics of the device contained 


within them, but these dynamics are not expressed on the system level. 
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1.2.1 Sparse Tableau Method 


The Sparse Tableau method is a very general method for describing systems 
employing the branch description of devices. Proposed in [4] and used in the ASTAP and 
SPICE simulators [1][15][16], the Sparse Tableau method breaks the system equations and 
variables each into three groups. The three sets of variables are the branch currents, branch 
voltages, and the nodal voltages. The three sets of equations are the node Kirchhoff 
Current Law (KCL) equations in terms of the branch currents, Branch Voltage equations 
relating branch voltages to nodal voltages, and the Constitutive equations relating branch 


voltages to branch currents. 


Figure 1.2.1-1 RC Example: Sparse Tableau 





Figure 1.2.1-1 shows an example of a simple RC charging circuit. Using the Sparse 


Tableau approach, the system variables are: 


i Voltage Source branch current 
es Resistor branch current 

Ic Capacitor branch current 

Vs Voltage Source branch voltage 
Pr Resistor branch voltage 

Vo Capacitor branch voltage 

e; Node | potential (voltage) 

e Node 2 potential (voltage) 


Note the reference node 0 is assigned a potential of 0. 
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The KCL equations are given by: 


The Branch Voltage equations are: 
Ve—e€,=0 Vo —€,=0 


Ve — (€,—€,) =0 


The Constitutive equations are: 
Woe yee : Ae 
Vr i ip = 0 
While the Sparse Tableau approach generates a consistent set of network equations, 
the size of the system is relatively large (eight equations in eight unknowns for this 


example). Furthermore, the method employs the branch description of devices which 


complicates the development of electro-mechanical models. 
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1.2.2 Modified Nodal Analysis 


The Modified Nodal Analysis method generates a compact set of system equations 
for systems of device models using branch descriptions. Modified Nodal Analysis (MNA) 
was formalized in [6], described in [1]{16], and employed in the circuit simulator MSINC. 
The procedure consists of writing the KCL equations for all but the reference node in terms 
of the branch currents, replacing the branch currents wherever possible with constitutive 
equations in terms of the branch voltages, appending any remaining constitutive equations, 


and substituting the branch voltages with the corresponding nodal voltages. 


Figure 1.2.2-1 RC Example: Modified Nodal Analysis 





Figure 1.2.2-1 shows a simple example of a simple RC charging circuit, the MDA 
approach would first write the KCL equations: 


Substituting the constitutive relations results in: 


a 1 ; 
lp +—v, =0 =—V,+¢ —-—0 
R Reve = alt 
The extra constitutive equation is given by: 


v,-V,=0 
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Substituting the nodal voltage results in the system of three equations 


adie ] de 
Aare: (ako -F(e,-€,)+C—-=0 


e,-V,=0 


While the Modified Nodal Analysis Method generates a compact set of equations, it 
does require the use of the branch description as well as the explicit definition of flow 


variables. Both restrictions can complicate teh modelling of electro-mechancial devices. 
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1.2.3 Standard Load Flow 


The Load Flow approach is traditionally used in the analysis of commercial power 
systems. For this application, the flow variables are usually real and reactive power while 
the potential variables are the voltage magnitude and phase angle. The Load Flow 
approach is a variation of nodal analysis described in many papers and texts on power 
systems including [14] [29] [31] [35] [49] [50] [76] [101]. The terminal description of 
devices is used since power flow is not conserved on the device level (The power entering 
a transmission line is not the same as the power leaving the same line). The basic 
procedure is to write the KCL equations in terms of the node potentials. Nodes with ideal 
potential sources are treated specially since their corresponding flow variable is not a 


function of the device voltages. 


Figure 1.2.3-1 RC Example : Load Flow 
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Figure 1.2.3-1 shows a simple RC charging circuit using the terminal description of 
the devices. A load flow approach using currents as the flow variable would result in the 


following procedure: 


Write the KCL Equation at nodes without potential sources 


lp tie, =9 


Substitute Constitutive relationships for the flow variables 


J A(Ve1 — Ve) 
p Re Ve) tC oe 0 


Ole 


a> 





Substiute the nodal potentials 


] de, 
a) eM 


All the remaining variables can be calculated from the solution of this differential 
equations. The load flow method definitely creates a very compact set of equations (only 
one in this case) but requires the flow variables be defined explicitly in terms of the 
potential variables, and must treat ideal potential sources as special exceptions. Neither of 


these restrictions is attractive for a general electro-mechanical system simulator. 
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1.2.4 WAVESIM Terminal Description 


The method employed in this thesis is similar to Modified Nodal Analysis with the 
exception that terminal potentials are used instead of branch voltages and that the 
constitutive equations are only expressed on the device level and never expressed on the 
system level. Potential difference equations are appendended to the system of KCL 
equations to equate explicitly defined potentials with their node potentials. For the RC 


example, the system variables are given by: 


is; Voltage Source terminal 1 current 
iC Capacitor terminal 1 current 

eo Node 0 potential (voltage 

e; Node | potential (voltage) 

a, Node 2 potential (voltage) 


Figure 1.2.4-1: RC Example: Terminal Description 





The KCL Equations for the RC example are given by: 
Isp + 8p ini (C1,€2) = 0 
lor + 8p inz(€1,€2) = 9 


Ig + 85 jso(lsy: Car 8c icollc; ,€)) =0 
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The Potential Difference Equations are given by: 
©. — &ociXicy» Co) = 9 
€1— 8s _ys1(ts;>€o) = 9 


€o- 8G lig) = 0 


Note that a reference device allowing for a more general method of setting the 
system reference points is employed rather than a reference node. While the number of 
equations is twice that of the Modified Nodal Analysis method, flows need not be 
conserved on the device level. Furthermore, the system of equations is easily partitioned 
into a sequence of five blocks for a more rapid solution (two 1x1 blocks, followed by a 
2x2 block, followed by two more 1x1 blocks). 
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1.3 Solving System Equations 


As stated earlier, the standard approach to simulating a physical system is to generate 


a system of differential algebraic equation of the form: 
Grit) 
O= g(x, y,u) 

To solve this system, it must first be converted to a system of purely algebraic 
equations by substituting the differential equations with discrete approximations. The time 
history of a variable is expressed as a series of discrete points in time where dynamics are 
expressed as algebraic relationships between the values of a variable at different discrete 


times. Standard methods for performing this approximation include the forward and 


backward Euler, Trapezoidal rule integration and Runge-Kutta methods. 


The major problem with this approach is the dependence of the time step on the fastest 
mode (smallest eigenvalue) of the dynamic system. This forces the entire system be solved 
with a very fine discretization of time, even though large portions of the system are not 
affected by the fast mode. 

In any case, the system of nonlinear algebraic equations must be solved. The two 
classes of solvers most commonly used are variations of the Newton-Raphson Method and 


several relaxation methods. 
1.3.1 Newton-Raphson Method 


The Newton-Raphson method works well for most systems as long as the initial 
guess for all of the variables are within the convergence region of the final solution. This 
method is used in SPICE and ASTAP and is based on a Taylor series expansion of the 


system of equations: 
F(x) =0=F(x,) +J(x,)A, +... 
See =X; +A, 


A,=—J '@,)FG,) 


The matrix J is called the Jacobian Matrix and its inverse must exist for the 


method to work. 
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1.3.2 Relaxation 


Relaxation methods assign one of the system variables to each of the system 
equations. After initial guesses are made for each of the variables, the variables are 
updated by solving their corresponding equation assuming none of the other variables have 
changed. The two most popular relaxation methods are the Gauss-Jacobi (popular with 
parallel processing computers) and the Gauss-Seidel method (usually used with serial 
processing computers). The Gauss-Jacobia calculates updates for all the system variables 


before actually performing the update: 
F(x) =0 
Ob aeeowe Teer eke cats eee = 0) 
The Gauss-Seidel method updates the system variables as the updates are calculated: 
F(x) =0 
i (ee camer Sree Ronee aeallp) = 0 
1.3.3 Waveform Relaxation 


An alternate method to solving the dynamic equations system wide is to solve them 
equation by equation over a given time interval. The Waveform Relaxation method 
represents variables by a sequence of points representing the time history of the waveform 
over a given time interval. Each variable can be discretized differently and is assigned one 
of the system equations. The system equations are solved over the waveform interval for 


their assigned variable with the other variables held at their current waveform values. 


Waveform Relaxation works well with loosely or directionally coupled systems, but 
does not work well for tightly coupled systems. The method does however, have good 
multirate performance since each differential equation can be solved using a time 


increment appropriate to it. 
1.3.4 WAVESIM Approach 


To summarize the traditional solving methods, the standard methods employing 
Netwon-Raphson can handle tightly coupled systems but perform poorly with multirate 


systems while waveform relaxation performs poorly with tightly coupled systems but 
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efficiently solves multirate problems. Unfortunately, the shipboard systems have both 
multirate and tightly coupled properties. For this reason, WAWVESIM combines the 


Newton-Raphson method with a waveform representation of variables. 


In WAVESIM variables are represented over a time interval by a vector of 
coefficients along with a type indicator for specifying how the coefficients should be 
interpreted. Common interpretations include Legendre Series coefficients, Chebyshev 
Series coefficients, and polynomial seres coefficients. For these representations, 
integration and differentiation are linear matrix operations and the issue of numerical 
stability of an integration technique disappears. Waveforms can usually be converted 


from one type to another with a linear matrix operation as well. 


With variables represented as vectors of coefficients, the Newton-Raphson method 
can be employed for solving tightly coupled systems. Good multirate performance 1s 
achieved through the linear matrix operator for integration along with waveform 


smoothing to average out phenomena faster than the time scale of interest. 
1.4 Thesis Outline 


This thesis focuses on developing a digital computer simulation environment suitable 
for studying shipboard electric power systems. WAVESIM, a simulation program written 
in the C programming language demonstrates algorithms for simulating systems of 
nonlinear lumped parameter models representing the electro-mechanical components 


composing an IED system. The key features of WAVESIM are: 


Devices defined independent of the encompassing systems 
Devices can be developed and tested without an exact knowledge of the 


topology of the systems incorporating the devices. 


Devices described using the Terminal Representation of devices 

Device constitutive relationships are written in terms of the actual values of the 
terminal potentials and not in terms of relative potentials. In this manner, device 
equations can be written in terms of a system reference when such a reference level is 
unambiguous. Furthermore, the flow variables are not required to be conserved on a 
device level. This greatly eases the task of modelling flows which also depend on a 


reference potential (power for example). 
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Devices defined independent of the manner in which terminal interface variables 
are expressed. 

Devices can be developed without specifying how the interface variables are 
specified. In WAVESIM, variables can be represented many different ways, all of 
which are irrelevant to the specification of the constitutive equations making up the 


device. 


System equations instead of the devices resolve input-output conflicts. 
WAVESIM does not constrain normal terminals where energy is transferred 


from having more than one output hooked together at a node. 


Interface Variables represented by waveforms 

Waveforms are a vector of coefficients which specify a given variable over a 
given time interval instead of a single value describing the variable at a given point in 
time. The waveform type determines how the coefficients should be interpreted for 
generating values of the variable within the time interval. Representing variables as 
waveforms has the primary benefit of removing the issue of numerical stability of 
integration techniques from the simulation. Integration and differentiation are merely 
operators on waveforms, no different from addition, subtraction, or any of the 


trigonometric operators. 


Differentiation and Integration performed on the device level instead of the 
system level. 

Most circuit simulators as described in the previous sections solve the 
differential equations associated with device constitutive equations on a system level. 
This method eases the task of evaluating the stability of linear systems but introduces 
new problems. If the eigenvalues of a dynamic system are widely separated in value, 
the simulation time step must be made very small for the entire system if conventional 
integration techniques are employed. WAVESIM solves the differential equations on 
the device level and employs waveform smoothing to remove dynamics which occur 


faster than the time scale of interest. 


While many of the pieces of WAVESIM are not new, several key concepts are 


presented in this thesis for the first time: 


The Terminal Description of devices 
Instead of specifying the interface of devices by ports consisting of a potential 


difference (branch voltage) and the flow through the potential difference (branch 
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current), the terminal description of a device assigns a potential and a flow entering 
the device for each normal terminal. Simulators based on branch voltages and 
currents require all of the flow entering a device to also leave the device. In this 
sense, the flow is conserved. The terminal description however, does not require 
conservation of flow within a device (Conservation of flow as expressed by 
Kirchhoff’s Current Law - KCL is required at connection points called nodes). The 
ability to construct models which do not conserve flows can simplify models where 
energy transformations occur, the reference potential 1s clearly known for the system 
and not just for the device, and certain forms of energy are not of interest. In many 
mechanical simulations for example, the amount of energy lost in friction is not of 
interest to the modeler. A simulation model based on branch potentials and flows of a 
device experiencing friction would be required to reject the frictional heat through one 
of its branches. 

The terminal description also allows for the transfer of information between 
devices through information nodes and information terminals. This feature is 
essential for successfully modelling many control algorithms. The ability to mix 
control signals and energy transfer through flow variables within the same simulation 


environment is a major advantage of the terminal description. 


The Structural Jacobian method for building and reducing systems 

The concept of the connection matrix for specifying the participation of system 
variables in system equations 1s expanded to include the structural form (1.e. diagonal, 
linear, nonlinear, etc.) of the dependence of the system equations on the system 
variables. The codes for the structural Jacobian adhere to a simple set of algebraic 
rules which can be used to construct a system structural Jacobian matrix from the 
individual device structural Jacobians. The system structural Jacobian facilitates the 
reduction of the numerical effort required to solve the system by identifying and 
characterizing a set of smaller blocks which when sequentially solved, determine all 
of the system variables. The system structural Jacobian can also be used to detect 


unconnected systems and indicate possible potential reference problems. 


The Systematic Treatment of Waveforms as an abstract data type 

WAVESIM departs from the conventional paradigin of representing variables in 
a dynamic simulation by a series of discrete points in time with a new paradigm based 
on representing variables as a sequence of waveform intervals. Within each 


waveform interval, the value of the waveform can be directly determined for any time 
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based on a vector of coefficients, a waveform type indicator for specifying how the 
vector of coefficients should be interpreted, and the time boundaries of the waveform 
interval. Devices are defined independent of the waveform type of the terminal 
variables. The principle advantage of using waveforms 1s that integration and 
differentiation are simple operators. The integral of a waveform is just another 
waveform. Simulation time steps are no longer controlled by the requirement for 
numerical stability of the integration technique. Instead, series truncation error 
contro] becomes the primary concem of the simulation environment. The ability to 
use arbitrary waveform types and convert between types allows the modeler to use the 


most appropriate waveform representation for the modeling problem. 


This thesis is composed of six chapters including this introduction. Chapter Two 
describes in some detail the specific properties of current shipboard electric power systems 
and proposed integrated electric drive systems. Chapter Three provides a framework of 
theory for developing the simulation environment WAVESIM and is broken into five 
subsections. The first subsection details the Terminal Description method for modelling 
devices. The second subsection demonstrates how to interconnect device models into 
systems, construct the system structural Jacobian, and generate a sequence of blocks for 
solving the system equations. The third subsection covers the treatment of waveforms as an 
abstract data type. Solving the system of equations employing waveforms is detailed in the 
fourth subsection. The fifth and final subsection of the third chapter covers modelling 
techniques and considerations not covered in previous sections. The actual WAVESIM 
implementation of the concepts developed in the third chapter are described in the fourth 
chapter. The fifth chapter presents results of several simulations conducted with 
WAVESIM. The final chapter provides an assessment of the work presented here as well as 


possible future developments. 


The appendices support the main chapters. Appendix A 1s a glossary of terms used 
through out this thesis. Appendix B details some possible problems with using continuation 
parameters. Appendices C and D are Load Flow examples of the terminal description 
method. Appendix E provides examples of waveform types and a number of operators for 
them. Appendix F presents a number of models useful for conducting shipboard power 


System simulations. Finally, Appendix G details the program files making up WAVESIM. 


This thesis introduces a number of new terms. To assist the reader, the first occurance 
of a new term is indicated by the distinctive Helvetica typeface. MATLAB variable names 


and sample sections of C programs are printed in Courier. 
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Chapter 2 Shipboard Electric Systems 
2.1 Typical Shipboard Electric Distribution System 


The electric power systems onboard naval warships differ considerably from the 
integrated power utilities found in developed countries. The differences arise from the 
small size of the shipboard systems and contrasting standards for optimization. Shipboard 
systems are optimized for survivability and minimization of weight and volume. Power 
utilities on the other hand optimize for reliability and minimization of cost. The unique 
characteristics of the shipboard systems result in markedly different design requirements 


and standards as compared to power utilities. 


Frigates, destroyers and cruisers are relatively small warships with corresponding 
small electric power systems. Frigates normally displace from 2000 to 4000 long tons 
(1 long ton = 2240 lbs) and have a primary mission of escorting merchant convoys. In the 
U.S. Navy, frigates have only one propulsion shaft and about half the armament of a 
destroyer. Destroyers displace from 4000 to 7500 long tons and are designed as escorts for 
aircraft carrier battle groups. Cruisers are larger than destroyers, displacing from 6000 to 
16000 long tons, carry more weapons, and are used to provide aircraft carrier battle groups 
with integrated anti-aircraft and anti-cruise missile defenses. U.S. Navy cruisers and 


destroyers all have two propulsion shafts. 


The installed electric plant capacity for U.S. warships has varied from 3000 KW to 
4500 KW per propulsion shaft over the past twenty years. Generally, the newer ships have 
more installed capacity. Figure 2.1-1 shows the electric plant characteristics for the major 
classes of conventionally fueled frigates, destroyers and cruisers constructed in the past 
twenty years. All the listed ships with the exception of the Knox class frigates use 
mechanically coupled gas turbine propulsion. The Knox class frigate is the last class of 
conventionally fueled warships to use 1200 psi steam for main propulsion. (All nuclear 
powered ships use 600 psi steam). Most of the Knox class fngates are presently being 


transferred to the reserve forces or being decommissioned. 
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| Ship Class (Nbr) Generator Type | Year 
FF-1052 Frigate Knox  [4x750KW| 3 Steam Turbine |1969| 
(46) 1 Diesel 
FFG-7 Frigate Oliver Hazard| 4x 1000 4 Diesel ] 
(1) Perry K 
DD-963 Destroyer Spruance 3 x 2000 3 Gas Turbine | 1975 
(31) K 
(4) K 


W 
WwW 

DDG-993 Destroyer 3x 2000 | 3GasTurbine [1981 
Ww 

DDG-51 Destroyer | Arleigh Burke | 3 x 3000 3 Gas Turbine | 199] 
(1 + 28) KW 


CG-47 Cruiser Ticonderoga 3 Gas Turbine 
(19 + 8) "Aegis" KW 





Figure 2.1-2 Shipboard Electric Distribution System 
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Figure 2.1-2 shows a typical ring bus architecture found on modern warships. The 


small size of the shipboard system results in many differences with respect to commercial 


systems. As a consequence the analysis of the shipboard plant requires recognition of these 


differences: 


l. 


ie 
8. 


Power Quality requirements relaxed relative to commercial 
standards. Constant frequency and voltage assumptions 


can not be made. See section 2.2 for more details. 


Very little Rotational Inertia require fast controls to 


maintain frequency. Infinite bus assumption does not hold. 


Transmission lines are very short and for the most studies, 


can be ignored. 


No scheduling of real or reactive power. All generators 
are loaded in equal proportion to their rating. 


Load Flow solution has little meaning. 
Load sharing information communicated to all online generators. 


Large loads (relative to the size of generation plant) present. 


Start up transients (load dynamics) are important. 
Power Electronic Switching loads are significant. 


Load shedding strategies are minimal. 


Figure 2.1-2 also indicates the requirement for a simulation environment to include 


the ability to model more than just electric power phenomena. Modelling shipboard 


systems also requires extensive representation of mechanical dynamics as well as energyless 


information transfer between components. This requirement is significant in that simulation 


packages for commercial power systems do not include this capability as an integral part of 


the simulation environment design. 


2.2 Shipboard Electric Plant Standards 


The primary standards for designing a shipboard electric plant are contained in the 


following references: 


Department of Defense, Interface Standard for Shipboard Systems, Section 300A, 
Electric Power, Alternating Current (Metric), MIL-STD-1399(NAVY), 13 
October 1987. 
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Department of the Navy, General Specifications for Ships of the United States Navy, 
Section 300, General Requirements for Electric Plant, Naval Sea Systems 
Command, 1987. 


Department of the Navy, General Specifications for Ships of the United States Navy, 
Section 320, General Requirements for Electric Power Distribution Systems, 
Naval Sea Systems Command, 1987 
The goal of electric power utilities 1s to provide a reliable source of high quality 

electric power at minimum cost. Shipboard systems on the other hand are designed to 

provide a survivable and continuous source of electricity. Quality and cost are secondary 


issues. Figure 2.2-]1 summarizes the minimum quality of power a shipboard system must 


provide 


Figure 2.2-1 clearly demonstrates the quality of power guaranteed onboard a warship 
is considerably lower than the quality of service provided by power utilities. Figure 2.2-1 
does not show however, how often the transient conditions occur. This information is 
provided by MIL-STD-1399 and summarized in figure 2.2-2. A major ramification of the 
low quality of power provided by the ship service electric system is that loads must be 
designed to operate and survive wide ranges of voltage and frequency fluctuations. This is 
one of the reasons why commercial equipment often can not be directly installed onboard 
ships (Shock requirements are also a major factor). Sensitive loads must provide their own 
filtering and protection circuitry. This militarization of equipment can add considerable 


cost and complexity to warship design, outfitting and maintenance. 
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| Figure ype Shipboard Electric Power Quality Standards (MIL-STD-1399) 
| ___ Frequency | 
Nominal 60 Hz 
Tolerance ama ey/o 
Modulation’ 0.5 % 
Transient Tolerance BI 0 


Transient Recover Time 2 seconds 
Worst Case Excursion £5.07 
Voltage 
Nominal 440/115 Volts 
Tolerance of 3 Phase Ave +5% 
Tolerance of any 1 Phase +7% 
Line Voltage Unbalance” 3% 
Voltage Modulation 2 Te 
Transient Tolerance + 16% 
Maximum Departure Voltage from +6% 
combination of 3 Phase Ave. and 
Voltage Modulation 
Worst Case Excursion + 20 % 
Recovery Time 2 Seconds 
Voltage Spike” 2500 / 1000 Volts 


| Voltage Waveform 


Maar Lone Deteition® 
Max Single Harmonic 
Max Deviation Factor . . | 
| Emergency | | | | 
Frequency Excursion -100 % to +12 % 
Voltage Excursion -100 % to +35 % 
Duration | =. | 2 Minutes P 





“100 measured over a period of 1 to 10 seconds. 
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1 Modulation (percent) = 


2 Line Voltage Unbalance is the difference between the largest line to line voltage and the 
smallest line to line voltage divided by the nominal voltage. 


3 A Voltage Spike is a voltage change of less than 1 ms duration. 


4 Total Harmonic distortion is the ratio of the rms value of the residue (after elimination of 
the fundamental) to the rms value of the fundamental. 


5 Deviation Factor is the ratio of the maximum difference between corresponding ordinates 
of the waveform and an equivalent sine wave to the magnitude of the equivalent sine wave. 
The equivalent sine wave is defined as having the same frequency and rms voltage as the 
wave being tested. 


535. 





Figure 2.2-2 Shipboard Electrical Reliability 


Voltage Transients of 10% or less Several times an hour 


Voltage Transients of 10% to 16% Several times a day 
Voltage Spikes above 200 Volts About once every 3 hours 





The basic reason for the low quality of power onboard ship is the lack of rotational 
inertia in the power system. In the commercial sector, the inertia of all the generators in the 
network add up to such a large number that no single fault can cause a frequency 
disturbance system wide. Onboard ship however, generators are often operated 
independently. Other than the inertia provided by motors, the only source of rotational 
inertia is the one generator. Since the generators are not very large, sudden load changes 
and faults can cause significant disturbances. Although speed governors and voltage 
regulators have improved significantly in the past twenty years, there is presently no way to 


prevent the transients from happening. 


The frequency tolerance limits in the steady state are rarely ever approached in 
modern warships. The rather loose tolerances allowed the use of droop governors to stably 
share loads. The electric plant operator on older ships could increase the load on a 
paralleled generator by increasing the base frequency set point on the mechanical speed 
governor. Adjusting the system frequency without changing the load sharing ratios required 
adjusting the base frequency set points on all the generator speed governors. On modem 
warships, all the generators normally operate isosynchronously and perform load sharing by 
transmitting load current information to Governor Control] Units which provide feedback to 


the isosynchronous governors. 
2.3 Shipboard Electric Plant Design 


In the commercial sector, the design of electric generation and transmission capacity 
are done continuously. Ships on the other hand, have a finite life (typically thirty years) and 
the expense of upgrading the capacity of the electric plant and distribution system once the 
ship is built is usually prohibitive. In this sense, capacity expansion onboard ships is not 


done. Instead, excess capacity 1s initially installed to account for projected growth in load. 


The maximum load for a ship design is determined by tabulating every load in an 
Electrical Load Summary and summing up the power requirements under different 
operating conditions. The maximum projected load usually occurs when the ship is in battle 
condition and the ambient temperature 1s low (Electric heaters are used in many areas of a 


ship). To account for uncertainty in estimating loads, a 20 % margin is added to the 
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maximum projected load. Another 20 % margin is added for capacity expansion 
requirements. Ninety percent of the capacity of all but one of the installed generators must 
meet or exceed the margined maximum projected load. The ninety percent requirement 
allows for imprecise load sharing when at maximum load while the all but one requirement 


accounts for taking one generator off line for maintenance. 
Figure 2.3-1 


U.S. Ships - Electrical Loads 
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Once the size of the electric plant is determined, there are a number of other 
considerations that must be accounted for. GENSPECS® require the system be ungrounded 
and based on Split Plant Operation (Each generator operating independently) with the 
capability for parallel operation. Electromagnetic Interference (EMI) and Electromagnetic 
Pulse (EMP) requirements place further constraints on the electric plant design and are 
detailed in MIL-STD-461 and MIL-STD-1310. Since warships are designed for combat, 
they must also be capable of surviving severe mechanical shocks from exploding ordnance. 
The shock requirements are particularly important for electrical equipment such as circuit 


breakers and generators. Specific requirements for shock are listed in MIL-STD-901. 


6 General Specifications for Ships of the Unites States Navy 





A number of loads onboard a ship are very important for survival of the ship and crew 
during combat and emergencies. These loads are designated vital loads and must be 
provided with primary and alternate sources of power. Some of the vital loads have 
automatic bus transfer switches (ABT) which switch to the alternate source automatically 
on loss of the primary source. Others use manual bus transfer switches (MBT). Examples 


of vital loads include: 


Collective Protection System Class W ABT 
Ventilation 
Emergency Communications MBT 
Emergency Lighting ABT 
Fire Pumps ABT 
AFFF Pumps ABT 
Interior Communications ABT 
Machinery Space Circle W Ventilation MBT 
Steering Gear Auxiliaries ABT 
Surface Search Radar MBT 
VHF Bridge-to-Bndge Radio MBT 
Vital Propulsion Auxiliaries MBT and ABT 
Auxiliaries to support generator prime MBT 
movers 


From a naval architectural viewpoint, the placement of electric generators requires a 
number of compromises. Placing the heavy generators as low as possible is beneficial for 
hydrostatic stability purposes. The lower the generator however, the more volume is 
required for intake and exhaust ducting. Gas turbine generators are lighter than diesel 
generators, but require greater volumes of air. Furthermore, design requirements exist for 
Separating 50 % of the installed capacity by two watertight bulkheads and installing a 
minimum of three generators. Generally, weight can be minimized by using the smallest 
number of generators (three). However, if four generators are used, the generators can be 
located in two machinery spaces instead of three. By using only one set of intake and 
exhaust ducts, volume for ductwork can also be reduced. Since most recent ships have had 
weight constraints placed on them by Congress, the minimum number of generators have 


been used.’ 


7 A very simple cost model for warships assigns a cost per ton of different components of a 
Ship. With this in mind Congress has in the past placed constraints on the weight of ships in 
order to keep costs down. 
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Enclaving is a concept for arranging ships which involves locating all the equipment 
required for a given combat system within the same general area of the ship. If a ship is 
completely divided into a number of enclaves, one enclave can be damaged by enemy 
ordnance while the others remain functional and capable of continuing the engagement. To 
work properly, this concept requires the enclaving of sources of distributed services (such as 
electricity, cooling water, fire fighting water and dry air). Presently, enclaving has not been 
incorporated in any warship design but its use has been proposed for a number of new 
designs®. If enclaves are ever used, they will have a significant impact on the type, size, 
number, and location of electric generators. In some enclaves it may not even be possible to 
locate a conventional generator. Alternate generating or storage devices such as fuel cells 


or batteries may be used. 


2.4 Integrated Electric Drive 


Most modem warships mechanically couple the main propulsion prime movers with 
the propeller shaft. The mechanical power train is very efficient but imposes constraints on 
machinery arrangement and adversely impacts survivability. The prime mover is usually 
very heavy and must be located near the center of the ship to prevent excess trim. Shafting 
must therefore penetrate a number of watertight boundaries and maintain precise alignment 
over a great distance. The long length and precision requirements of the shafting make it 
very vulnerable to weapon induced damage. While electric propulsion eliminates many of 
the survivability and arrangement constraints of the mechanical system, the propulsion 
system must be carefully designed to ensure overall plant efficiency is not degraded by the 
extra power Conversion losses in converting to and from electric power. Designed properly, 
an electric drive system can achieve the survivability and arrangeability benefits without 


suffering from a lower propulsion plant efficiency. 


Integrated electric drive interconnects the generation of power for propulsion with the 
generation of ship service electric power. The propulsion plant for U.S. warships typically 
averages between 30 and 37.5 MW per shaft. The capacity is sized to provide enough 
power to propel the ship at a desired maximum speed. Most ships however, do not operate 
for extended periods of time at maximum speed. Operating at half maximum speed requires 


only about 20 percent of the installed power and quarter maximum speed requires only 2 or 


8 Enclaving requires a greater redundancy of equipment which results in the ship becoming 
larger and more expensive. Since most ship designs are cost constrained, enclaving 
provisions are often deleted to reduce the per unit price of the warships. 
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3 percent. Thus a 28 knot fngate with a 30 MW plant could go 7 knots using less than 1 
MW of power and 14 knots with about 6 MW of power. If the propulsion plant consists of 
two 15 MW generators, one generator could easily supply all the required power for both 
propulsion and ship service at the normal operating speeds of 12 to 15 knots. This has the 
potential of reducing the fuel consumption of warships under normal operating conditions 
by improving the overall efficiency of the power plant even though the efficiency of the 
power transmission system is lower. By careful selection of generator number and size, one 
can tune the overall efficiency of a plant for optimization at several different speeds. In the 
U.S. Navy, optimizing plant efficiency for 20 knots is beneficial since this is the speed used 


to calculate the amount of fuel carried by the ship.’ 


In a typical integrated electric drive scheme, the propulsion prime movers are 
connected to both a propulsion generator and to a ship service generator (PDSS or 
Propulsion Derived Ship Service). The speed of the generator is set to optimize efficiency 
of the prime mover at the given power loading. Consequently, cycloconverters are used to 
convert the power to either 60 Hz for ship service, or to whatever frequency the propulsion 
motors require. Usually, an additional diesel or gas turbine ship service generator 1s 
included to provide power in port or during emergencies. Figure 2.4-1 shows a typical 


PDSS design for a two shaft frigate sized ship. 


Figure 2.4-1 emphasizes the need to model mechanical dynamics and control 
information signals. The control signals can couple the dynamics of different devices 
within the system and must therefore be carefully modelled. The control signals can also 
destroy such properties as diagonal dominance which makes analysis of commercial power 


systems much easier. 


One of the features of an electric drive system which may be exploited in the future is 
the ability to divert all of the propulsion power capacity from propulsion to some sort of 
high power combat system. Weapons such as rail guns and high energy lasers may become 
possible. These types of weapons would be safer for the ship since the requirement to store 
large amounts of chemical explosives for propellent charges would be reduced. Energy to 


move projectiles would be stored in the form of relatively inert fuel oil instead of highly 


9 Most other navies use 18 knots which allows for combined plants such as CODOG where a 
diesel engine is used for cruising and a gas turbine for high speed. Unfortunately, the size 
requirement for a diesel capable of propelling a ship at 20 knots is prohibitive and results in 
U.S. warships only using gas turbines and carrying much more fuel. 
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Figure 2.4-1 Integrated Electric Drive 
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explosive chemical propellents. Switching large amounts of electric power onboard ships 
presents a number of technical challenges both in the design of physical equipment and also 
in attempts to accurately simulate the phenomena. The effect of pulse loads on the electric 


system 1S not a trivial simulation problem. 
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Chapter 3 Framework 


Conducting time domain simulations of systems of nonlinear lumped parameter models 
characterizing shipboard electric power systems requires an organized approach to 
developing device models as well as network equations. The major contribution of this 


thesis is the development of a simulation environment having the following properties: 


1. An object oriented approach to developing the mathematical description of devices 


independent of the manner in which the variables are represented. 


2. An organized method for generating system equations for interconnecting device 


models into subsystems and systems. 


3. An algorithm for solving the system equations and variables by identifying smaller 
blocks of equations and variables which can be sequentially solved. The 
algorithm develops the concept of the device structural jacobian matrix and the 


system structural jacobian matrix. 


4. The ability to use a wide range of methods to descnbe variable waveforms. In 
particular, describing waveforms through vectors of coefficients of polynomial 


series, orthogonal function series, and data series are stressed. 


5. The ability to solve the system of equations by employing either the 
Newton-Raphson Method or Waveform Relaxation. The Newton-Raphson 
method is modified to improve convergence properties through the use of 


continuation methods. 


This chapter is organized into five parts. The first part defines the device which is the 
fundamental building block of the system simulation. The second part shows how to 
interconnect several device models into systems and subsystems. The third part defines the 
waveform as a vector of coefficients to approximate waveforms over time intervals. The 
fourth part details the actual procedure for conducting a simulation. The fifth and final part 


details some finer points which should be considered when constructing models. 
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3.1 Device Description 


A Device Description is an organized manner for describing the charactenstics of 
a physical component. This description includes definitions of variables which interface 
with other components in a system, variables called states which allow for information 


storage, and constitutive relations descnbing the device behavior. 
3.1.1 Interface Variables 


The interface variables are defined as either potential variables or flow 
variables depending on their interaction with the interface variables of other devices 
within a system or subsystem. Systems and subsystems are constructed by grouping 
the interface variables of one or more devices into sets called nodes and applying 


network equations determined by the types of variables attached to the nodes. 


All potential variables attached to a node are equated to a potential value 
associated with the node. Physical quantities which can be classified as potentials include 
voltages, signal levels, rotational speeds, deflections, and pressures. All potentials are 
referenced to 9. All potential variables connected to the same node must be defined with 
respect to the same system wide reference level. In other words, 0 must mean the same 


thing for all of the potentials attached to a given node. 


The sum of all flow variables attached to a node is equated to zero. Physical 
quantities analogous to flow variables include currents, power flows, torques, forces and 


mass flow rates. 
3.1.2 Terminals 


Terminals provide a mechanism for organizing the interface variables of a device. 
In general, there are two types of terminals: Normal Terminals and Information 
Terminals. 


A normal terminal has associated with it a flow variable and a potential variable. 
Its electrical analog is one of the wiring terminals on an electrical device. A mechanical 
analog is the rotating shaft coupling of a gearbox. The equations for exchanging energy 
between devices can be generated through the list of normal terminals connected together 


at a given node. 





An information terminal has associated with it only a potential variable. The 
potential variable is used to convey knowledge between devices without transferring 
energy. Set points, meter readings, and control signals are all examples of energyless data 


which can be conveyed through information terminals. 


All normal terminals have an associated KCL Group number. A KCL group is the 
smallest subset of a device’s terminals such that the sum of the flow variables within the 
subset is identically zero for at least one of the possible dynamic configurations of the 
device. Normal terminals which can not be associated with a KCL group are given a group 
number of 0. The remaining terminals are assigned the group number of their parent KCL 
group. 

The KCL Group number Is used to detect possible reference frame problems within a 
simulation network. A given electrical circuit problem for example, must have at least one 
normal terminal with a 0 group node within a given independent system to ensure the set 
of system KCL equations is not singular. Normally this terminal is associated with a one 
terminal device with an export potential and import flow which is used to specify the value 
of a given reference node potential. This Reference Frame Check is discussed in greater 


detail in section 3.2.4. 


Some devices may have variable numbers of KCL Groups depending on the 
operating point of the device. A simple model of a two terminal switch for example, 
would have 1 KCL group when the switch is closed (the sum of the currents entering the 
switch is identically zero) and 2 KCL groups when the switch is open (both flow variables 
are identically zero). For the purpose of defining the device, the worst case in terms of 
creating singular systems should be used. In the switch example, each terminal should 


have their own KCL group number for a total of two KCL groups. 
3.1.3 Variable Direction: Import and Export Variables 


The Interface variables can further be classified by whether they are a resource 
(Import) or product (Export) of the device description. A device description can be 
considered a means for generating export variables based on the values of the import 


variables, states, parameters, continuation parameter, and time. 


An import variable is taken as input by the device description. An import variable 


can be any interface variable associated with either normal or information terminals. To 
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ensure a consistent set of equations when several devices are connected together in a 
system, the total number of import variables associated with normal terminals must equal 


the number of normal terminals 


An export variable is explicitly defined and considered a product of the device 
description. An export variable can be any interface variable associated with either normal 
or information terminals. To ensure a consistent set of equations when several devices are 
connected together in a system, the total number of export variables associated with 


normal terminals must equal the number of normal terminals. 
3.1.4 States 


States are variables whose values are stored for a given time for later use. States 
can be used for example, to store the constant of integration for a dynamic equation. States 
can also be used to store the operating mode for a given device. In general, if the value of 
a given variable depends on the previous value of another variable, that other variable is a 


state. 
3.1.5 Parameters 


Parameters are constants which specify characteristics of the device or in other 
words, customizes a given device description to represent a given physical device. A 
model of a resistor for example, includes a parameter for resistance. This precludes the 
requirement to develop a model for every resistor value. We only need construct a generic 


resistor model instead of a 10K resistor model, a 22K resistor model, etc. 
3.1.6 Constitutive Equations 


The constitutive equations are a consistent set of equations for specifying the 
values of the states and export variables. In general, the number of constitutive equations 
needed is equal to the number of normal terminals plus the number of export variables 
associated with information terminals. The total number of import variables associated 
with normal terminals and the total number of export variables associated with normal 
terminals must independently equal the number of normal terminals. There is no 


constraint on the number of import variables associated with information terminals. 
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3.1.7 Device Jacobian Matrices 


A Device Jacobian Matrix provides the sensitivities (partial derivatives) of the 
export variables with respect to the import variables. This implies there is a given 


ordering of both the import X, and export X, variables: 
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The Device Jacobian Matrix is used to generate a consistent set of import variables 
which simultaneously satisfy the device constitutive equations along with constraints 
imposed by the connections of terminals to nodes. From the device point of view however, 


the Jacobian matrix 1s merely a product that must be computed. 


Up to this point, we have not discussed the manner in which the variables are 
described. If the variables are represented by real numbers, then each element of the 
Jacobian is also a real number. If instead the variables are represented by vectors, then the 


Jacobian elements will be matrices. 


3.1.8 Device Structural Jacobian Matrix 


The Device Structural Jacobian Matrix describes the properties of the elements 
of the device Jacobian matrix for a given type of variable representation without actually 
providing any values. The following codes can be used to describe the properties of the 


matrix elements of the device Jacobian matrix: 


Type of Matrix 


Zero Matnx (all e elements are 2 eae eae 
Identity Matnx (always the identity matrix) 
Diagonal Matrix (always a linear main diagonal matrix) 


Linear Matnx (The elements are always constant) 


Nonlinear AC Matrix (see Note 3.1.8-1) 


Nonlinear Matrix (The elements may not be constants) 


Unknown (The dependence is unknown (treat as nonlinear)) 





Note 3.1.8-1: An AC Matrix is one for which the constant component of the export 
variable depends only on the constant component of the import variable. The other 
components of the export variable can not depend on the constant component of the 
import variable but are not restricted in any other way. 

The device structural Jacobian matrix is useful in developing the algorithm for 
generating a consistent set of import variables without having to deal directly with the 


potentially much larger device Jacobian matrices. If an iterative solution scheme is used to 
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develop the consistent set of import variables, the device structural Jacobian matrix 
indicates directly which matrix elements must be recalculated for each iteration. (Only the 


nonlinear and unknown elements have values which change between iterations) 
3.1.9 Continuation Parameter 


A system containing one or more nonlinear devices may be difficult to solve with an 
iterative method. The region of convergence around the solution may be so small as to 
make the probability of success for choosing a starting point for the iterative scheme 
almost zero. One method for enlarging the region of convergence is through the use of a 
continuation parameter which varies from 0 to 1. When the continuation parameter 
has value 1, the export variables are developed using the normal nonlinear constitutive 
equations. When the continuation parameter has value 0 however, the export variables are 
developed using a linear set of constitutive equations. As the continuation parameter 
increases from 9 to 1, the export variables traverse a continuous path from the linear 
solution to the nonlinear solution. One common method for generating such a dependence 


on a continuation parameter © Is: 
F(X, Q@) = OF (X)+(1-Q)F(X) 


where F(X) is the nonlinear function for generating the export variables, F{X) is the 
linear function approximation, and F(X,q@) is the function for determining the export 
variables for intermediate values of @ Section 3.4.2 descnbes in detail continuation 


parameters in relation to the Newton-Raphson method. 
3.1.10 Discontinuity Time Prediction 


If the variables are described as a waveform over a given time interval [f,.f,] 
knowledge of the time of discontinuities can prove useful to the algorithm which generates 
the consistent set of import variables. The accuracy of a vector descnption of a waveform 
often deteriorates greatly if there is a discontinuity during the time interval. Varying f, 
such that it falls on a discontinuity will often improve the accuracy of the waveform 
representation. For this reason, each device has the opportunity to recommend a 
recalculation time for the current interval. Normally, the system would use the minimum 
recommended recalculation time offered by any of the devices to recompute the time 


interval. 
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3.2 Network Description 


A network is composed of a system of devices and subsystems whose terminals 
are interconnected at nodes. The network is a closed system having no terminals defined 
for any of its nodes. A Subsystem is a system having terminals defined for at least one of 


its nodes and therefore can not be solved independently of other devices or subsystems. 
3.2.1 Nodes 


A node connects together one or more terminals from one or more devices. The 
nodal connections are the means by which devices are combined to form systems (both 
networks and subsystems). The nodes provide the association of device import and export 
variables with system variables through nodal equations. Each node is assigned a 
serial number for identifying it from the other nodes. There are two types of Nodes: 
Normal Nodes and Information Nodes. 


3.2.1.1 Normal Nodes 


A Normal Node has at least one normal terminal attached to it. Information 
terminals can be associated with the node as long as none of the information terminal 
potentials are defined as an export variable. A normal node has associated with it a node 
potential as well as a Kirchhoff Current Law (KCL) equation. The number of normal 
nodes is designated by n,. 


In a subsystem, a normal node can also have associated with it a terminal for 
connecting with other subsystems and devices. This terminal can be either a normal 
terminal having an associated terminal potential and flow variable or an information 
terminal having only an export potential. (import and export refer here to the direction 
relative to the defining subsystem which 1s opposite to the normal definition which is 
relative to the components of the subsystem). The total number of normal node normal 
terminals defined for a subsystem is designated n,,,. For any given subystem the number 
of normal node terminal export variables and import variables must both independently 


equal n,,,- The total number of normal node information terminals is designated n.,,,. 
3.2.1.2 Information Nodes 


An Information Node has only information terminals attached to it. Furthermore, 
one and only one of the terminal potential variables must be an export variable. Only a 


node potential is associated with an information node. Information nodes work in the 
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same manner as hooking up stereo componenents: you can hook up as many inputs 
(import variables) as you want to any given output (export variable), but should never 
hook up two or more outputs together. The number of information nodes is designated by 


N,. 


As an option for subsystems, an information node can have associated with it an 
information terminal for connecting with other subsystems. Since the meanings of import 
and export are once again reversed for this terminal, no other export potentials from other 
devices or subsystems may be attached to the node if the information terminal potential is 
an import variable. If the information terminal potential is an export variable, exactly one 
other export potential from other devices or subsystems may be attached to the node. The 


total number of information node information terminals 1s designated 7,,,. 
3.2.2 System Variables 


System variables comprise the minimum set of variables from which all of the 
device import and export variables can be derived from. The set of system variables 1s 
composed of node potentials as well as all device import flow variables and normal node 
normal terminal export flow variables. For a subsystem, the node terminal import 
variables are assumed to be provided by the encompassing system or subsystem and are 


not considered system variables. 
3.2.2.1 Node Potentials 


All of the node potentials of the normal and information nodes are system variables 


which must be solved for. Hence there are a total of n, = n, + n; node potentials. 
3.2.2.2 System Flow Variables 


All of the Import Flow Variables of the various devices making up the system as 
well as the export flow variables of the normal node terminals are system variables. The 


number of system flow variables is designated by ny. 


3.2.3 System Equations 
3.2.3.1 Kirchhoff Current Law Equations 


Kirchhoff’s current law states the sum of the flow variables entering a node is equal 


to zero. For a given normal node or normal terminal node, this law is expressed by 
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generating a list of the terminals of the various devices and subsystems attached to the 


node. The number of Kirchhoff Current Law equations is equal to the number of normal 


nodes n,. 
FO = x 1, =0 
where 
Sf KCL Equation for node j (Should Equal Zero) 
n, Number of normal terminals attached to node 
I; Flow Variable associated with ith normal terminal attached to node J 


3.2.3.2 Potential Difference Equations 


A Potential Difference Equation is created for each of the export potential 
variables of the various devices and other subsystems and for each of the import potential 
variables of the node terminals. This equation merely states the difference between the 
node potential and the potential variable 1s zero. This equation is expressed by 
generating a list of the terminals of the various devices and subsystems attached to the 
node having an export potential variable. Since one and only one export information 
potential can be assigned to an information node and can never be attached to a normal 
node, the number of potential equations due to export information potentials is simply 7;. 
The requirement for a device to have equal number of import and export variables 
associated with normal terminals forces the number of export normal potentials to be my. 


Hence the total number of potential equations 1s n, = n; + ny. 


FiO =Vj-Vi=0 


where 
FiO Potential Difference Equation for node j export potential variable 1 
(Should Equal Zero) 
V; Node j Potential 
Vi; ith export potential variable associated with node j. 
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5.2.3.3 R win and Ge 


One method for preventing linear dependences among the system equations is to 
modify the equations to include an extra term corresponding to either a small conductance 
Gain to the ground potential for KCL equations or a small series resistance R,,,,, for the 


potential difference equations. The KCL equation is now given by: 


f0=G,,V,+ X1,=0 


i=] 
The potential difference equation is similarly modified: 


FiO =Vj—Vie—Ronnd g = 0 


min j} 


The goal in using G,,;, and R,,;, 1S to reduce the condition number of the system 
Jacobian matrix to the point where the system can reliably be solved (A singular 
matrix has an infinite condition number). G,,,, and R,,;, can also add fictitious dynamics 
to the system and thereby lead the simulation to produce incorrect results. Hence if used, 
G,,in and R,,;, Should be large enough to bring the condition number down to a reasonable 
level, but small enough to prevent their inclusion from having appreciable effect on the 


simulation results. 
In general, the use of G,,,, and R,,;, Should be avoided for these reasons: 


l. G,,;, and R,,;, are fictitious elements. If either is significant, they should be 


explicitly included as a device. 


2. The indiscrimuinite use of G,,,, and R,,;, adds to the complexity of the system 
and decreases the degree to which the system can be reduced into smaller blocks. 
In other words the inclusion of G,,,, and R,,;, may greatly increase the computation 
time. 

Gin and R,,;, are included in WAVESIM for these reasons 


1. Gain and R nin Can be selectively specified for individual nodes. If a simulation 
fails to converge for one reason or another, G,,;,, and R,,,, can be employed to find 
the part of the system experiencing difficulties. G,,, and R,,;,, are excellent 


debugging tools. 


Soe 





2. Since G,,,, effectively connects the node to the ground potential, G,,,, can be 
used to ensure all of the nodes have the same potential reference and ensure there 


are no linear dependent KCL equations. 
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3.2.4 Reference Frame Testing 


If a given set of a system’s normal nodes can be found such that all terminals 
attached to any of its nodes have nonzero KCL groups and such that if a terminal is 
attached to one of the set’s nodes, then all of remaining terminals of the parent KCL group 
are also attached to one of the nodes of the set, then there exists the possibility of a singular 


system due to the linear dependence of the KCL equations for the set of normal nodes. 


If G,,, 18 non-zero for a node, it should be considered a terminal with a 0 KCL 


Group. IfG,,,;, 1S zero, it should be ignored. 


Testing for a possible singular system can be accomplished with the following 
algorithm: 


1. Set all the normal node circuit group indicators to 0. 
Set the circuit group counter to 0 


Set the circuit group singular flag to 0 


2. Start with the first normal node having a0 circuit group indicator 
If none can be found then algorithm is complete. 


Increment circuit group counter. 


3. Change the circuit group_indicator of the node to the 


circuit _group_counter. 


4. For each terminal attached to the node: 

4a. If the KCL group number is zero, set the 
circuit group singular flagto i; 

4b. If the KCL group number is nonzero, loop through each normal 
terminal of the device. If the terminal belongs to the same 
KCL group and the node the terminal is attached to has a 
Q circuit group indicator, then set the node 
circuit group indicator to the negative 


of the circuit group_counter. 


5. Search all of the nodes for anegative circuit group indicator 
If none can be found and the circuit _group_ singular flag 1S zero 
Wam user that a singular system may exist with the group nodes. 
If none can be found then go to step 2 


If one is found, then go to step 3 
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Setting a proper reference for each such set of system nodes can be accomplished by 


attaching to one of the nodes a one terminal device having the following characteristics: 


3.2.4.1 Reference Device 


Interface Variables 


Terminal Potential Variable Flow Variable (KCL Grp) Type 
Ref V (export) I (amport) (0) Normal 
Parameters 
Veg Reference Potential Level 
Equations 

V= Vref 


Device Structural Jacobian 


Js = [0] 


Device Jacobian 


Jy = [0] 


Notes 


Most conventional circuit simulations define a reference node for which a potential 
is defined and the KCL equation is not written. Adding this reference device to a node 
effectively converts that node to a reference node in the usual senses. While it is true that 
the KCL equation and an additional Potential Difference equation are still written for this 
reference node, each is part of a one element block. The potential difference equation can 
be solved before the simulation starts since it does not depend on any of the system 
variables. The flow variable on the other hand, only appears in the KCL equation of the 
one node and thus can be solved after all the other system variables have been found. In 
fact, the flow variable should normally equal zero if the rest of the circuit is indeed 


linearly dependent. 


As a convenience to the user, WAVESIM automatically attaches a reference device 


with V,.,= 0 to the node with serial number 0 if that node is used. 
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3.2.5 System Reduction 


The previous sections detail a method for generating a full set of system variables 
and system equations. The total number of system variable equals n, =n, + n; + n, which 
also equals the number of system equations. For even a small system the algebraic order 
n, Can become quite large. For this reason, elminating system variables and equations 
through system reduction is desirable. The primary tool for performing system reduction 
is the System structural Jacobian. 


3.2.5.1 System Structural Jacobian 


The System Structural Jacobian facilitates the reduction of the algebraic order 
of the system by showing the nature of the dependence of system equations to each of the 
system variables. The System Structural Jacobian is constructed by combining elements 
of the device structural Jacobian matrices according to the arithmetic of structural 


Jacobian elements. The types of elements in the system structural Jacobian is given by: 


Nonlinear Matrix (The elements may not be constants) 


ee ee en: ~~ 


Unknown (The dependence is unknown (treat as nonlinear)) 





Note 3.2.5.1-1: An AC Matrix is one for which the constant component of the export 
variable depends only on the constant component of the import variable. The other 
components of the export variable can not depend on the constant component of the 


import variable but are not restricted in any other way. 


The addition and subtraction operators for the structural Jacobian elements is a 
function of the manner in which the system variables are represented. For all of the 


methods used in this thesis, the following definitions apply: 
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I+l=D 
I-IT=0 
’ Mame | Ge f 
-I+0=D 
-[-[=D 
tntm=tmine=n(n2m, nz) 
U>N>A>L>D>I>90 


Note, the Identity Code I, is not strictly necessary and if eliminated simplifies the 


addition and subtraction operators to: 
tnhtm=tmtn=n(n2m) 


Before the system structural jacobain can be constructed, the system variables and 
equations must be ordered. The first 7, variables are the node potentials of the normal 
and information nodes arranged in the order of the node serial numbers. The next n, 
variables are the import flow variables ordered first by device then by device terminal. 
The first 2, equations conform to the Kirchhoff Current Law equations for the normal 
nodes arranged in order of the node serial numbers. The remaining n, equations are the 
potential equations for the export potentials ordered first by the node serial number they 
are attached to, then by the order of the devices attached to the node, and finally by the 


order of the terminals in the device. 


The system structural Jacobian is constructed in two parts after being initialized to 
contain only 0. First, a Kirchhoff Current Law equation is generated for each normal 
node. The normal terminals of the normal nodes are examined one at a time. If the flow 
variable is an import variable, it is also a system variable and an J is added to the 
corresponding element of the system Jacobian matrix. If the flow variable is an export 
variable, its corresponding row of the device structural Jacobian matrix is extracted. The 
columns of the device structural matrix row correspond to the device import variables. 
All of the device import variables can be associated to either a node potential (one of the 
first nz, columns of the system structural Jacobian) or to one of the remaining n, import 
flow variable columns. Hence it is quite easy to locate to which column each element of 
the device structural Jacobian row must be added. If G,,, is non-zero for the node, a D 


code is added to the column corresponding to the node potential. 
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The remaining n, rows of the system structural Jacobian matrix are constructed by 
examining each node one at a time. If the node has an export potential associated with it. 
An I is added to the corresponding node potential column and potential equation row 
element (unless of course the node is a reference node and does not have a column 
associated with its potential). The row of the device structural Jacobian matrix 
corresponding to the export potential is then extracted. In exactly the same manner as 
described above for the export flow variables, the columns of the system structural 
Jacobian matrix are correlated to the columns of the device structural Jacobian matrix. 
Once correlated, the elements of the device structural Jacobian row are subtracted from 
the appropriate elements of the system structural Jacobian matnx. If R,,,, 1s non-zero for 
the node and the terminal having the export potential has an import flow variable, then a 
D is added to the column corresponding to the import flow flow variable. If R,,;, 1S 
non-zero for the node and the terminal having the export potential has an export flow 
variable, then a D is multiplied by the elements of the corresponding row of the device 
structural Jacobian matrix before being added to the corresponding column in the system 


structural Jacobian matrix. 


Once the structural Jacobian matrix has been constructed it can be examined to 
ensure there are no glaring problems such as a row or column containing only 0 elements. 
If a row or column contains only 0 elements, the system is ill-posed and can not be 


solved. 
3.2.5.2 Blocks 


The primary reason for constructing the system structural Jacobian matrix 1S to 
break down the system of equations and system variables into smaller blocks which can 
be sequentially solved instead of solving the entire system at once. A block B; is defined 
as n,; System variables and n,; equations which only depend on system variables of the 
present block and previous blocks in the sequence. A block of size n,, 1s identified by 
finding n,; rows in the system structural Jacobian matrix that have not already been 
allocated to a block and have exactly n,; columns containing non-0 elements. Of the 
many combinations of blocks which can be found for a system, the best combination 


contains the largest number of small blocks. Here is an algorithm for finding the blocks: 


1. Create a list for each row containing the number of unallocated non-0 entries 


in that row. (Initially all the rows and columns are unallocated) 
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2. Examine the list for rows having only | unallocated non-0 entries. Create a 
block for each of these rows and their associated columns. Mark the rows and 


columns as allocated. 
3. | Update the list of unallocated non-0 entries in each row. 
4. Continue steps 2 and 3 until no more single rows can be allocated. 


5. Examine the list for two rows only having unallocated non-@ entries in the 
same two columns. Create a block for each pair of rows and their associated 


columns. Mark the rows and columns as allocated. 
6. Update the list of unallocated non-@ entries in each row. 


7. Repeat steps 2-6 until no more single row and double row blocks can be 
identified. 


8. Examine the list for three rows only having unallocated non-0 entries in the 
same three columns. Create a block for each set of three rows and their 


associated columns. Mark the rows and columns as allocated. 
9. Update the list of unallocated non-0 entries in each row. 
10. Repeat steps 2-9 until no more blocks of up to size 3 can be identified. 


11. Continue the above algorithm until all of the rows and columns have been 
allocated. Remember it is necessary to go back and attempt to identify 
smaller sized blocks after discovernng a larger block since the removal of a 


column could allow the identification of anew smaller block. 


The order of identifying blocks is very important because they must be solved in the 
same order. Each block contains the same number of system variables and system 
equations. The equations only depend on system variables determined from the present 
and previous blocks. Hence the simulation problem becomes an issue of solving 


sequences of relatively small systems of equations described by blocks. 
3.2.6 Reduced System 


The reduced system consists of the sequence of blocks which when solved, provide 
the solution for all the system variables. Solving each of the blocks can be done a number 
of ways. Most schemes start with an initial guess for the system variables and generate 


corrections to the guesses until all of the system equations for that block are satisfied. 
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Generating the corrections is normally done through the use of a block Jacobian matrix 
which can be constructed in much the same manner as the system structural Jacobian. If 
the block structural Jacobian does not contain any A, N or U elements, the block Jacobian 
can be inverted and multiplied by the system equation errors to provide the required 
corrections. If there are any nonlinearities, this scheme can be performed several times 
until the system equation errors are close to zero. This method is commonly referred to as 
the Newton-Raphson method and if the initial guess is close enough to the solution, the 
method converges quadratically. This method is described in much more detail in section 
3.4.1. 


Relaxation techniques can also be used to calculate the system variables. Relaxation 
techniques start with an initial guess for all of the system variables and update each 
variable one at a time by solving a single system equation by assuming all of the other 
variables are constant. Typically, one system equation 1s assigned the task of solving for a 
particular system variable. With careful thought as to the assignment of variables to 
equations, it is often possible for such a system to converge to a solution. Common 


relaxation techniques are the Gauss-Seidel and Gauss-Jacobi methods. 
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3.3 Waveforms 


Up to this point, the development of the simulation structure has been independent of 
the manner in which variables are actually described. The simplest and most commonly 
used method for representing variables is through a single real number representing the 
value of a variable at a specific tume. For static simulations where the problem is to obtain 
the steady state solution for the system, this method works very well. Appendix C and 
Appendix D demonstrate this procedure for the classic load flow problems. For dynamic 
simulations however, some knowledge as to the time history of the variables is needed to 
calculate derivatives and integrals. A dynamic simulation is implemented as a series of 
static simulations where the dynamics are represented by functions of the time increment 
and state variables. The various integration techniques for this type of simulation differ 
only in the interpolation scheme used to approximate the variables between successive static 
simulations. The time increment between static simulations must be carefully controlled to 
ensure the interpolation scheme has enough accuracy for numerical stability. Integration in 
this manner requires careful control of the time increment to ensure the interpolation 


scheme is accurate enough to ensure numerical stability along with an accurate solution. 


Another approach to representing variables is the waveform. This method employs a 
vector of coefficients to continuously describe the time domain value of the variable over 
some time interval [4 ¢,]. The type of the waveform determines how the coefficients are 
interpreted to generate the time domain values. Possible types include Data Series, Fourier 
Series, Legendre Series, Polynomial Series and Legendre Series. The principal advantages 


of using waveforms over discrete points include: 


iz Interpolation is not generally required to determine intermediate points. The 


value of a variable can readably be determined for any time. 


2. The numerical stability of Integration and Differentiation techniques do not 
have to depend on the time step control since integration and differentiation 
become waveform operators on an equal level to all other operators. Time step 
control becomes only an issue of numerical accuracy and not of numerical 


stability. 


Se Certain operations may be easier to perform with one waveform type. The 
ability to efficiently convert a waveform from one type to another type and 
back again allows one to use the most efficient waveform type in the 


calculations of a given operator. 
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3.3.1 Waveform Definition 


A waveform approximates the instantaneous value of a variable over some time 


interval. The elements of information contained within a waveform must as a minimum 


include: 


4. 
>. 


The name of the waveform 


The beginning and ending times of the interval (f, f,) 


An Array of Coefficients representing the waveform (c;) 


The number of coefficients in the Coefficient Array (7) 


A waveform type indicator. 


The waveform type indicator identifies how the coefficients should be 


interpreted when operations are performed on the waveform. Here is an example of a C 


structure defining a Waveform: 


typedef struct Waveform 


{ 


} 


char *name; 
double EO; 
double ales 
void ai el- 
long n; 
Long Eype, 


long version; 


struct Waveform 


hel 


/* 
/* 
/* 
/* 
{* 
/* 


character string of the name 
of the variable */ 
time of the beginning of the interval */ 
time of the end of the interval */ 
array of coefficients */ 
number of elements in the array */ 
waveform type indicator */ 
Version of this waveform */ 


*next; /* pointer for forward 


lanked ivste */7 


Struct Waveform *last; /* pointer for backwards 


linked lists */ 


SEruct Jacobian *jnum: /* pointer to Jinked list of 


jacobians where this waveform 
is the numerator */ 


Struct Jacobian. *j;den; /* pointer to linked list of 


WAVEFORM; 


jJacobians where this waveform 
is the denominator */ 


The above definition also includes the following optional information: 


mee ee a gen 


A Version Number to record a change in the waveform’s properties. 


An Address Pointer to the waveform representing the previous time interval. 


An Address Pointer to the waveform representing the following time interval. 


An Address Pointer to a linked list of Jacobian Stnuictures. 
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The waveform address pointers allow one to construct a linked list of waveforms to 
describe the time history of a variable over a number of time intervals. The Jacobian 
structure as well as the version number will be described in section 3.3.3. 


Note the waveform coefficients are declared to be of type void. This is done to allow 
for the coefficients to be abstract data representations in themselves. Normally the 
waveform coefficients would be double precision floating point numbers, but it should also 
be possible to incorporate other types of data. It may be advantageous for example, to 
represent the coefficients with complex numbers. In this case, each element in the 
coefficient array would be a structure holding double precision floating point numbers 
corresponding to the real and imaginary parts (Or magnitude and phase angle) of the 


complex number. 





3.3.2 Waveform Operators 


Waveform Operators are functions which act on waveform arguments to generate 
new waveforms, or provide some information about the waveform arguments. The types 


of functions can be broken down into several groups: 


I. Arithmetic Operators 


eZ: Tngonometric/Exponential Operators 
a Switching Operators 

4. Integral/Differential Operators 

a: Waveform Content 


6. Special Functions 
3.3.2.1 Arithmetic Operators 


The arithmetic operators are the customary addition, subtraction, multiplication, 
division, and assignment operators usually associated with floating point arithmetic. The 
assignment operator is a bit more complex since it must incorporate waveform type and 


number of coefficient conversions. 
3.3.2.2 Trigonometric/Exponential Operators 


The Trigonometric/Exponential operators include most of the transcendental 
functions used in engineering. Examples include sine, cosine, tangent, logarithms, 
exponentials, as well as the inverse functions. Error handling can become quite complex 
since several of these operators may be undefined at one or more points within the 
argument waveform. These operators are usually handled by converting the arguments to 
a series of data points, performing the operation point by point, and then converting back 


to the appropnate waveform type. 
3.3.2.3 Switching Operators 


Switching Operators are operators producing waveforms which themselves or one 
of their derivatives are discontinuous. Examples include the absolute value function, the 


sign function and the step function. The typical method for calculating these functions is 
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to determine the discontinuity points and use integration to create a characteristic function 
series solution (e.g. Legendre Series or Chebyshev Series) for the result. The series 


solution is then converted to the appropriate waveform type. 
3.3.2.4 Integral/Differential Operators 


One of the key advantages of using waveforms in dynamic simulations is that 
integration and differentiation become very simple operators where the stability of a 
numerical integration scheme 1s generally not an issue. For many waveform types, the 
integration operator is a linear matrix operation with bounded coefficients. If the 
argument waveform has bounded coefficients, the retumed waveform will also be 
bounded. Of course, numerical stability does not assure numerical accuracy. Because the 
integration operator typically generates some truncation error, the returned waveform can 


stil] contain considerable errors. 
3.3.2.5 Waveform Content 


The significance of the Truncation Error of a waveform can be estimated by 
calculating the waveform content of its higher order term. The waveform content of a 
term is defined as the magnitude of a coefficient divided by the square root of the sum of 
the squares of all the coefficients. Normally, one expects the higher order terms of an 
orthogonal series representation to progressively have smaller and smaller waveform 
contents. Hence if the last few terms have values below a preset threshold, the truncation 


error can normally be assumed negligible. 


Accurate truncation error estimation 1s still a difficult and currently unsolved 
research topic. The waveform content method is a practical method but should not be 


taken as the last word on the subject. 
3.3.2.6 Special Operators 


Several special operators unique to waveforms should also be developed. One very 
useful operator returmms the time of zero crossing of the waveform. Another returms the 


value and time of every local minimum and maximum of a waveform. 


The smoothing operator is one method for reducing the waveform content of 
higher order coefficients. A waveform is smoothed by returning the local average of the 


waveform over some prespecified time increment. Smoothing eliminates discontinuities 
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in a waveform and its derivatives. Since discontinuities tend to amplify the waveform 
content of the higher order terms, removing the discontinuities should reduce the higher 


order term waveform content. 
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3.3.3 Jacobians 


A Jacobian matrix contains the partial derivatives of the coefficients of one 
waveform with respect to another waveform. Here is a sample C structure to define a 


Jacobian: 


typedef struct Jacobian 
{ 

struct Waveform *num; /* address of waveform in the 
numerator of the partial 
derivatives */ 

struct Waveform *den; /* address of waveform in the 
denominator of the partial 
derivatives */ 


long version; /* Version number of the 
jacobian matrix */ 
long num version; /* Version nbr of 
numerator Waveform */ 
Jong den version, j/* Version nor of 
denominator Waveform */ 
Vvord.%* 3; /* array of jacobian elements 


The first row index is for 
an array of pointers whose 
elements are arrays with 
the colum index */ 

ehar Ss); /* Structural Jacobian Code */ 


struct Jacobian *next; /* address for linked list 
of Jacobians */ 
} 


JACOBIAN; 


Jacobians are used in the process of solving simultaneous systems of waveform 
equations through relaxation methods or through the Newton-Raphson Method. The 
purpose of num version and den_version is to record which versions of the numerator 
and denominator waveforms the jacobian was calculated for. The element version is used 
when several jacobians are combined and it is necessary to determine whether the 


combined matrix must be recalculated. 


In general, all operations defined for a waveform should also generate the jacobian of 
the results with respect to the arguments. Through the use of the chain rule, the jacobian 
matrix of the export variables of a device with respect to the device import variables can be 


determined. 


The structural jacobian code indicates the dependence and structure of the 


jacobian matrix. Here is a list of the codes: 
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Type of Matrix 


Zero Matrix (all elements are always zero) 


Identity Matrix (always the identity matrix) 


Diagonal Matrix (always a linear main diagonal matrix) 
Linear Matnx (The elements are always constant) 
Nonlinear Matnx (The elements may not be constants) 


Unknown (The dependence is unknown (treat as nonlinear)) 





The structural Jacobian code along with the version numbers determines whether or 
not a jacobian matrix needs to be recalcuated. If the structural jacobian is of type 0, I, D, 
or L then the jacobian need not be reconstructed if the there is a version mismatch between 
the waveform version and the jacobian version. If the structural jacobian of type N or U, 
and there is a mismatch between the version numbers of the jacobian and the waveforms, 
then the jacobian elements must be recalculated. After every recalculation, the version 
numbers are updated. In this manner, only Jacobian matrices with changing coefficients 


are ever recalculated. 


Technically, the structural jacobian codes depend on the waveform type used. In this 


thesis however, all of the waveform types produce the same structural jacobian codes. 
3.3.3.1 Jacobian Operators 


Several operators for jacobian objects will prove useful in developing a simulation 


environment. These operators include: 


Addition and Subtraction 

Identity and Zero Jacobian generators 
Multiplication by a constant 
Multiplication of two jacobians 


Multiplication of a jacobian by a waveform 


Oy. Se age te se 


Inverting a jacobian 


If the waveform is described by an array of double precision floating point numbers, 
the Jacobian coefficients can also be defined to be an array of double precision floating 


point numbers. In this case, the above operations employ standard matrix manipulations. 
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3.3.4 Waveform Examples 


While the possibilities of waveform definitions is endless, this thesis will concentrate 


on the following waveform types: 


a 
oe 
3 









Legendre Series 


Polynomials 







Matlab Polynomials 


Chebyshev Series Tae 


The code in the above table refers to the value of element type in the WAVEFORM 
structure. Appendix E describes these waveforms and their arithmetic in great detail. 
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3.4 Conducting the Simulation 


Once the physical system has been specified by device descriptions and network 
equations, the solution for all of the system variables can be determined in several ways. 
The method used in this thesis is the Newton-Raphson method with continuation 


parameters. 


3.4.1 Basic Newton-Raphson Algorithm 


The Newton-Raphson method solves a system of nonlinear equations F(x,u) = 0, 
FQ e€ &’, for the system variables x € KR” with system input variables u € KR” by first 
linearizing the system of equations about a given guess for the solution x“ then solving the 
linear system to produce a new guess x**’. This procedure is repeated until F(x*,u) = 0 is 
satisfied within a given tolerance. The sequence of points x“ starting with k = 0 is called 
the solution trajectory for x°. A converging solution trajectory eventually converges to a 


solution while a diverging solution trajectory does not. 
F(x,u) is linearized by taking the Taylor series expansion about the point x*: 
F(x,u) = F(x",u) +J(x*,u)x,+O(x 1) =U 
ae to 
where the Jacobian matrix J(x,,u) 1s defined by: 


_ OF (x, u) 


UiGee u) oe 


Assuming the error O(xex) 1s negligible and the Jacobian can be inverted, the 


correction x, for a given guess x“ is given by the linear approximation: 
-l,_k k 
x,=—J ux 


The correction is applied to x* to produce x**", the value of x for the next iteration: 


x'tlex* +x, 
Around each solution of F(x,u) = 0 for which the Newton-Raphson method reliably 
converges, a region exists such that if a trajectory enters that region, it will never leave and 
eventually converge to the solution. The size of this local convergence region depends on 


the nonlinearity of the system. For purely linear systems, this region encompasses the 
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entire §” space. If the intitial guess falls within the local convergence region, the 
Newton-Raphson method will by definition converge. If the initial guess falls outside the 
local convergence region, one of several things can happen. First, the solution trajectory 
could enter the local convergence region of a solution and converge on a Solution. Second, 
the Newton-Raphson method could fail due to a singular Jacobian. Third, the trajectory 
could diverge and tend to infinity. Fourth, the trajectory could become cyclic where 
x**4 = x“ for k sufficiently large enough. Finally, the trajectory could enter a chaotic region 


in which there is no solution but from which the trajectory never leaves and is not cyclic. 


As an example, define F(x,u) to be the following 1X1 system: 
F(x,u)=x°-x 


Figure 3.4.1-1: F(x,u) = x°-x 
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The Jacobian matrix 1s: 


: re ae : 
The recursion formula for x°”’ is given by: 
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k+l _ Jk (x") —x" 
eae ara 
3(x*) - 1 
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There are three solutions for this system and their local convergence regions are 


given by: 


Local Convergence Region 
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In two other regions, the solution trajectory jumps to one of the local convergence 


regions after one iteration: 
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Convergence Region 





In two other regions, the solution trajectory may jump to one of the local 


convergence regions after several iterations or fail to converge: 


Variable Behavior Region 
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On the boundaries for the above regions, the Newton-Raphson method fails: 
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In the above eee no constraints were made in the speed of convergence or on 


2 
the size of x. If |.x* |» 1 the speed of convergence will be very slow since x**' = 5x" and 


the number of iterations / required will be about: 


_log(|x* |) 


= = 5.68] : 
log(1.5) 68 log(| x |) 


Furthermore, most machines have a limit as to the largest number which can be 
represented. If an iteration causes x to exceed this number in magnitude, a floating point 
overflow error will typically be generated. This phenomena is known as Newton Overflow 
and has the effect of reducing the size of the convergence regions. For example, if x is 
known to be bounded by the interval [-10 10], then x’ should be restricted to the following 


— -Converzence e Region 
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regions: 
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3.4.2 Continuation Methods with Newton-Raphson 


The previous discussion indicates the need for careful selection of the initial guess x”. 
The use of a continuation parameter in so called homotopy methods is one of the many 
ways for attempting to generate x° within the convergence region of the desired solution. 
In general, a function H(x,u,a)=90 is generated such that H(x,u,1) = F(x,u) and 
H(x,u,0) = G(x,u) where G(x,u) is a linear function in x. One common method of creating 


A(x,u,Q) 1s: 
H(x,u,a)=oOF (x,u)+(1-a)G(,u) 


The problem now is to develop the linear function G(x,uw). There are several 


approaches which can be taken for each row G,(x,u): 


1. Linearize about a known operating point. This is equivelent to providing an 
initial guess for each of the variables and using the Newton-Raphson method 
directly. 


2. Use a least squares fit of a linear system over a known operating region of 
FAx,u). 


3. Select Gf{x,u) such that the solution for H(x,u,0) = 0 is most likely to be within 


the convergence region of F(x,u). 
Once H(x,u,a) has been constructed, it can be used in several ways: 


1. Start with o=0 and obtain a solution to the linear system, then progressively 
increment alpha by small amounts and solve the nonlinear system until a=I. The 
rational is to employ the unbounded local region of convergence of the linear system 
to move the initial guess into the local region of convergence for the next nonlinear 
system formed by incrementing a. As @ is incremented, the solution for the previous 
value of a is assumed to be within the local region of convergence for the present 
value of alpha. Appendix B demonstrates this may not always happen due to 


bifurcations of solutions as @ 1S incremented. 


2. Start with a=1 and attempt to obtain a solution to the nonlinear solution. If the 
trajectory has not converged after n,,,, iterations, decrement and attempt to find a 
solution. Progressively decrement & until a solution is obtained, then increment @ 
using the solution of the previous value for a for the initial guess. This procedure 


assumes the local convergence region for a given solution will increase as q@ is 
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decremented. Eventually the local convergence region should grow large enough to 
encompass even a poor guess for the solution. This procedure has the advantage over 
the previous method in that it may avoid bifurcations which occur between 0 and the 


minimum value for @ used. However, the number of iterations for @ may be larger. 


Note that the value for ,,,, aS well as the convergence criteria may be a function of 
a. There is no reason to obtain a highly accurate solution for intermediate values of @ 
since the only purpose is to move the initial guess for the next & iteration into the new 
local region of convergence. Only when a=] should the convergence criteria be enforced 


for obtaining a highly accurate solution. 
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3.4.3 Simulation Algorithm 


The simulation algorithm employed by WAVESIM is conducted totally within the 
MATLAB environment and is composed of four parts. The first part initializes all of the 
simulation parameters. The second part performs the time increment control and has 
embedded with in it the third part which is the sequential solving of each of the blocks. 


The final part is composed mostly of plotting and storing the results of the simulation. 


Figure 3.4.3-1: Simulation Flowchart 
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3.4.3.1 System Initialization 


A number of parameters and arrays need initialization before the simulation can 


commence. These parameters and arrays are: 


sb n min 


sb n max 
sb _n dat a 


sb_dt_ 1 nit 


sb dt optimum 


sb dt min 


sb dt _max 


sb_dt_ ave 


Initial number of waveform coefficients 
Actual number of waveform coefficients used 


Waveform type indicator 


Beginning time of simulation 
Ending time of simulation 


Minimum number of coefficients to use 
Maximum number of coefficients to use 


Number of points per waveform for plots 


Initial time increment 
Optimum time increment 


Minimum time increment 
Maximum time increment 


Minimum time of interest (Averaging interval) 


Break Points are user specified times for which waveform interval boundaries are 


forced to occur. Break Points are completely optional and their inclusion is up to the 


system modeler. 


sb_ bp 
sb_bp_nbr 


sys node serial 


sys node name 


sb_alpha_init 


sb dalpha init 
sb dalpha_min 
sb dalpha_max 
sys Gmin 

sys Rmin 


Array of Break Points 
Number of break points 


Array of Node Serial Numbers 
Array of Node Names 


Initial Value of continuation parameter alpha for nonlinear 
blocks 
Initial Value of alpha increment 


Minimum alpha increment 
Maximum alpha increment 


Array of Gmin values for all of the nodes 
Array of Rmin values for all of the nodes 


The index for sys Gmin and sys_Rmin are the node numbers of the nodes they 


apply to. 
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sb_check_ eqn _err 


sb check var err 


sys_kcl_err 


sys pot_err 


sys _nd_err 


sys fv err 


sb_i_kcl_err 


sb_i pot_err 


sb ind err 


sb i fv err 


= Q for don’t check equation error 

= | for checking equation error 

= 0 for don’t check max variable correction 
= 1 for checking max variable correction 


Array of maximum KCL errors for all nodes 
Array of maximum Potential Differences for all nodes 


Array of max corrections to Node Potentials for all nodes 
Array of max corrections to Flow Variables for all nodes 


Multiplier for maximum KCL error 
for alpha less than 1 

Multiplier for maximum Potential Difference 
for alpha less than 1 


Multipher for max correction to node potential 
for alpha less than 1 

Multipher for max correction to flow variable 
for alpha less than 1 


The index for the above eight arrays are the node numbers of the nodes they apply 


to. 


sb_maxcnt 


sb _i_maxcnt 


sb_div_start_cnt 


sb_div_max_cnt 


sb i div_err 
sb_max wc 
sb_nbr we 
sb_mult_we 


sys pot _scale 


sys_flow_scale 


Maximum number of iterations for alpha = 1 
Maximum number of iterations for alpha < 1 


Number of iterations to skip before checking 
for divergence 

Maximum number of diverging iterations before 
assume syStem is diverging 

Multiplier of errors for ignoring diverging check 


Maximum waveform content of a waveform 
Number of coefficients to apply waveform content to 
Multiplier to sb_max we for decrementing N 


Array of Scaling factors for node potentials 
Array of Scaling factors for flows attached to nodes 


The index for sys pot scale and sys flow _scaie are the node numbers of the 


nodes they apply to. 


dev par_name 


dev sO name 


Device parameter arrays: name is the device name 
Device State initial value array: 
name is the device name 
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Ivar na nbs 


ivart nd _nbr 


ivar fv name 


ivart fv name 


his t 


his N 


his col 


his nd_nbr 


his fv_name 


his s_name 


blk_nbr_nrow 
blk nbr_ncol 
blk_nbr_row_sys 
blk _nbr_col_sys 


blk _nbr_linear_ flag 


Initial guesses for node potentials: 
nbr is the node serial number 
Waveform type for initial guess 
nbr is the node serial number 
Initial guess for flow variables: 
name 1s the variable name 
Waveform type for initial guess 
name is the variable name 


Matrix of time increment end points 
First row is beginning of intervals 
Second row is end of intervals 
Columns are waveform interval index 


Vector of number of coefficients in waveforms for each 
waveform interval 


The waveform interval index. After simulation this equals 
the number of columns in history arrays 


Matrix of Node Potential waveforms. Each column 
corresponds to the waveform for the node potential over a 
given waveform interval. nbr is the node serial number 


Matrix of Import Flow Variable waveforms. Each column 
corresponds to the waveform for the import flow variable 
over a given waveform interval. name is the variable name 


Matnx of Device name state values. The first column 
corresponds to the initial state values with subsequent 
columns corresponding to the state values at the end of 
waveform intervals. Note this matrix has 1 more column 
than all the other history arrays. 


Number of rows in block nbr 

Number of columns in block nbr 

Cross Reference of Block nbr rows to System Rows 

Cross Reference of Block nbr columns to System Columns 


= 0 if block nbr is nonlinear 
= 1] if block nbz is linear 
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Time Increment Initialization 


ddt Actual time increment 
tt0 Beginning of current waveform interval 
tei End of current waveform interval 


ddt, tt0, and ttl are initialized according to the following equations: 


ddt = sb dt init 
ttO = t0 


ttl = minimum of: 


tO + ddt 
tL 
sb _bp(1l) 
ent tot Set to zero: Total number of Jacobian inverses 
his flops Number of floating point operations used 
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3.4.3.2 Time Loop 


Truncation Error Control! 


The simulation time interval between t0 and t1 may be divided into a number of 
waveform intervals to improve the truncation error of the system variable waveforms. In 
general, truncation error can be reduced by either increasing n or by decreasing the 
waveform interval tt1 - tt0. Within WAVESIM, the general strategy for dealing with 
too large of a truncation error is to increase the number of coefficients n if the waveform 
interval is less than sb_dt_optimum and shorten the waveform time interval if greater 
than sb_dt_ optimum. In general, the strategy is to minimimize N while maximizing the 
waveform interval subject to the constraint that the truncation error is within tolerances. 
Finding the optimum combination of waveform intervals and number of coefficients 1s 


not obvious and much work remains for developing better algorithms. 


3.4.3.2.1 Time Loop iteration initialization 


The simulation time loop continues as long as tt0 < t1. The beginning of each 


iteration begins with the definition of the following arrays: 


tt = [ttO ttl sb dt ave] 
ii = Identity Matrix of size nN 
zz = Zero Matnx of size NxNn 


Variable Initial Guesses 


Next, initial guesses are provided for all system variables (var_nd_nbr and 
var fv_name) by converting the waveforms ivar_nd_nbr of type ivart_nd_nbr and 
waveforms ivar_fv_name of type ivart_fv_name into waveforms of type wtype and 


size N. 


In the present incarnation of WAVESIM, the same waveform is used as the initial 
guess for all waveform time intervals regardless of the values for ttO and ttl. 
Normally, a constant value is specified. A better method would allow the user to specify 
an actual guess as to the waveform history as a function of time. The time loop iteration 
initialization would then have the responsibility of converting the waveform data as 
provided by the user into a waveform of type wtype and size N over the interval between 
ttO and tt1. Providing an initial guess for the waveform history of all the variables 
would allow for example, a linear model of a system be run first to generate the initial 


guess for a nonlinear model of the same system. Convergence of the nonlinear system 
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should be greatly accelerated for many systems. Parameter sensitivity studies would 
also be greatly accelerated if the parameter variations are not expected to cause major 


changes in system performance. 


Failure Flags 


Two final variables, converge failure and fatal_error are initialized to zero. 
converge failure 1s set to one by a block if convergence failed for that block or if one 
of the block waveforms has too large of a harmonic content. Convergence could fail if 
the number of iterations exceeded the maximum allowed and the alpha increment is 
smaller than the minimum allowed. converge failure is used to indicate the following 
blocks should not be solved because previous blocks could not be solved. fatal_error 
is set to one 1f convergence cannot be obtained even when N is equal to or greater than 
the maximum value sb_n_ max and the time increment is equal to or smaller than the 


minimum value sb_dt min. If fatal error is set, the simulation fails. 
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3.4.3.2.2 Solving the Blocks 


The blocks are solved sequentially in the order of their detection in the system 
reduction procedure. If converge failure is nonzero, a previous block could not be 
solved for the given time increment and number of coefficients. For this reason, a block 


is not solved if converge failure is nonzero. 


Figure 3.4.3-2: Solving the Block 
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3.4.3.2.2.1 Block Initialization 


Each block requires the initialization of several arrays and variables before the 
block can be solved: 


84 - 





blk nbr_max_eqnerr 


blk nbr max varcor 


blk _nbr_ imax eqnerr 


blk nbr_imax varcor 


blk _nbric nt 


blk nbr_cnt_ div 


blk _nbr_ alpha 


blk nbr dalpha 


good_alpha 


good var_nd_nbr 


good var fv_name 


blk nbrt rec 


blk nbr_ive 


div cnt 


Giv err 


Array of maximum errors for the block equations 
Array of maximum variable corrections for the block 
variables 


Array of multipliers to blk_nbr_max_eqnerr for 
alpha < 1 
Array of multipliers to blk _nbr_max varcor for 
alpha < 1 


Number of iterations (initialized to 0) 
Number of diverging iterations (initialized to 0) 


Block continuation parameter. 

= ] if linear block 

= sb alpha_init if nonlinear block 

Block continuation parameter increment 

= sb dalpha init 

Last value of alpha for which block converged. 


Initialized to -1 


Last value of node nbr potential for which block 
converged. Initialized to var_nd_nbr 


Last value of import flow name for which block 
converged. Initialized to var_fv_name 


Recommended recalculation time for block 
Initialized to tt1 


Array of indexes in block variable array for which the 
variable was greater than allowed. 
Initialized to an empty array. 


correction 
Number of diverging iterations, set to 0 


Maximum relative error of previous iteration 
Initially set to 0. 
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3.4.3.2.2.2 Continuation Parameter Loop 


The block continuation parameter loop continues as long as blk_nbr_ alpha < 1. 


Within this loop, the following procedures occur: 


1. Import Variables for all associated devices specified 
2. Device Objects called to generate 
A. Export Variables 
B. Device Jacobian Matnx 
C. State values at time tt1 
D. Recommended recalculation time 
3. Block recalculation time calculated 
4. KCL and Potential Difference Equation Errors calculated 
5. Errors Scaled and compared to maximum limits 
if good, solution saved and blk_nbr_alpha incremented 
as necessary. 
6. Iterations counted and compared to maximum limit 
blk_nbr_alpha decremented and variables reset 
as necessary. 
7. Block Jacobian Matnx assembled and scaled 
8. Variable Corrections Calculated 
9. | System variables corrected 
3.4.3.2.2.2.1 Device Import Variable specification 


The matrix dev_i_name is generated for each device name where the columns 


are the waveform coefficients for each of the device import variables. Each column of 


the dev_i name matrix is one of the system variables, hence all are available. 


3.4.3.2.2.2.2 Call Device Objects 


Each of the device objects associated with the block 1s provided with the 


following information: 


wtype Waveform type 


dev_1 name Device name unport variable matrix 
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dev_par_ name Device name parameter array 
dev_sO_name Device name State initial value tt0 array 


tt Time structure 
blk_nbr_ alpha Block abr Continuation Parameter 


From this information, each of the device objects generates the following 
dev_e name Device name export variable matrix 

dev_j_ name Device name jacobian matrix 

dev_sl_name Device name State final value tt1 array 


dev _tr_name Device name recommended recalculation structure 
= [ntl ntt]) where 
ntl = recommended tt1 for present interval 
or Set to tt 1 1f no recommendation 
ntt = recommended tt1 for next interval 
or set to tt.0 if no recommendation 


3.4.3.2.2.2.3 Recommended Recalculation Time 


The block recommended recalculation time blk _nbr trec is Set to the 
minimum value of all the nt1 values from all of the devices associated with the block. 


If convergence fails blk_nbr_trec 1S used to generate a new value for tt1. 


Similarly, blk_nbr_ntrec is set to the minimum value of all the ntt values 
greater than tt1 from all of the devices associated with the block. For a successful 
convergence, blk_nbr_ntrec is used to help generate a new value for tt1 for the 


next waveform interval. 
3.4.3.2.2.2.4 Equation Errors 


For each of the node nd KCL equations associated with block abr, an error 
variable blk_nbr_kcl_nd is generated by adding the flow variables of the attached 
terminals to the flow through Gmin. Likewise, for each of the export potential name 
Potential Difference equations associated with block abr, an error variable 
blk _nbr_pot_name is generated by subtracting from the node potential waveform, the 
waveform of the export potential as well as the contribution from Rmin 


blk nbr_kcl_nd = x dev_e name(:,col) + var _fv_vwname + 
var_nd_nd X Gmin 
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blk _nbr pot _vname = var _ndnd - dev_e_name(:,col) - 


dev_x_name(:,col) xX Rmin 


where 
nbr Block Number 
nd Node Serial Number 
name Device name 
vname Variable name 


(:,col) | The appropriate column from the matrix 
x Either e or i depending on associated flow variable 
being an export or import variable 


The KCL equation errors are multiplied by the appropriate flow variable scaling factor 
from the sys flow scale array while the Potential Difference equation errors are 
multiplied by the appropriate potential scaling factor from the sys pot scale array. 


Once scaled, the error vectors are assembled into a block error vector blk_nbr_err. 


3.4.3.2.2.2.5 Error Criteria Check 
Applying Error Criteria 


If blk_nbr_alpha 2 1 then blk_nbr_ ier is filled with the indexes of the rows 
of blk_nbr_err which are greater in magnitude than the corresponding rows of 
blk nbr max eqnerr. In the same manner, blk_nbr_rel_err is set equal to the 


absolute value of blk_nbr_err divided by blk_nbr_max_eqnerr. 


If blk _nbr_alpha < 1 then blk_nbr ier is filled with the indexes of the rows 
of blk_nbr_err which are greater in magnitude than the corresponding rows of 
blk nbr_imax eqnerr. Similarly, blk_nbr_rel_err is set equal to the absolute 
value of blk nbr err divided by blk_nbr_imax eqnerr. 


Divergence Check 


On the first iteration for a given value blk_nbr_alpha, div_cnt is initialized to 
0. For the first sb div _start_cnt - 1 iterations, div_err is set to the maximum 
value of blk_nbr_rel_err. On subsequent iterations, if the maximum value of 
blk _nbr_rel_ err is smaller than div_err then div_cnt is reset to 0, otherwise 
div_cnt is incremented. In any case div_err 1s set to the maximum value of 


bikiobr rel err. If div_cnt 2 sb _div_max_cnt then the algorithm assumes the 
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block is diverging for the given value of blk _nbr alpha. The failure to converge 
condition 1s indicated by setting blk nbr_cnt = maxcnt: either sb i _maxcnt if 


blk_nbr_alpha < 1 0r sb_maxcnt if blk nbr alpha > 1. 


Block Convergence Success 


If blk_nbrier is the empty set or sb check eqn err is 0, and 
blk_nbr_alpha 2 1 and blk_nbr ive is the empty set, then the block solving 
algorithm has been completed and the continuation parameter loop is broken. The 
algorithm proceeds to checking the truncation error for the system variables associated 
with the block. 


Increment Continuation Parameter 


If blk_nbr_ier is the empty set or sb check_eqnerr is 9, and 
blk_nbr_alpha < 1 and blk_nbr_ivc 1s the empty set, then it is time to increment 
the continuation parameter blk_nbr_alpha. First however, the current value of all 
the variables associated with the block are copied into good var _nd_nd or 
good_var fv_name. blk_nbr_alpha is copied into good_alpha. The variables and 
continuation parameter must be saved because it may be necessary to restore the 
variables if the block fails to converge with the next continuation parameter value. 
blk_nbr_alpha is then set equal to the minimum of 1 and 


blk nbr_ alpha + blk_nbr_dalpha and the continuation parameter loop is repeated. 


Iteration Count: Decrement Continuation Parameter 


If the error is still too large, corrections to the system variables associated with 
the block must be generated. But first, the number of iterations blk_nbr_cnt must be 
incremented and compared to the maximum allowed maxent: either sb_i _maxcnt if 
blk_nbr_alpha < 1 OF sb_maxcnt if blk _nbr_ alpha 2 1. If the limit has been 
exceeded, and one of the devices has recommended a value for blk_nbr _trec less 
than tt1, then converge failure is set to 1 and attempts to solve the block cease. If 
the limit has been exceeded and blk_nbr_trec equals tt1, the block is recalculated 
with a decremented blk_nbr_alpha which is set to the maximum of: 

(blk_nbr alpha + good_alpha) / 2 
blk nbr alpha - blk _nbr_dalpha 
0) 
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If blk_nbr_alpha has been decremented, the system variables associated with 
the block must be reset to the values stored in either good_var_nd_nd or 


good var fv_name. 
Block Convergence Failure 


If blk_nbr_alpha - good_alpha < sb dalpha_min then the block has failed 
to converge and nothing more can be done on the block level. The variable trec is 
set equal to blk_nbr_trec and the converge failure flag is set to 1. This 1s a 
signal to the system to not solve any more blocks and either adjust the value of tt1 or 


adjust the number of coefficients n before trying to solve the system again. 


3.4.3.2.2.2.6 Assemble Jacobian 
Jacobian Construction 


If the error is too large, but the maximum number of iterations maxcnt has not 
been exceeded, the block Jacobian matrix must be calculated. The block jacobian 
matrix blk_nbr_3 1S constructed in the same manner as the system structural jacobian 
was previously constructed with the exception that now the variables and equations 
are only those which are part of the block and the matrix elements are submatrices 


instead of structural jacobian codes. 


Jacobian Scaling 


Once the block jacobian has been assembled, it 1s scaled by dividing each of the 
columns by the appropriate element of either the sys_flow_scale (if the column 
corresponds to an import flow variable) or sys_pot_scale (if the column corresponds 
to a node potential) vectors. Likewise, rows of the block jacobian are multiplied by 
the appropriate element of either the sys_flow_scale (if the row corresponds to a 
KCL equation) or sys_pot_scale (if the row corresponds to a Potential Difference 
equation) vectors. Scaling is performed to normalize all of the variables and 
hopefully improve the accuracy of the numerical computations required for solving 


the variable corrections. 


Correction Vector Calculation 


The variable correction vector blk_nbr_dlta 1s generated by solving the matrix 


equation: 


90 - 





blk_nbr_j blk_nbr dita = blk_nbr err 


The most direct method (and one of the least numerically efficient method) of 
calculating blk_nbr_dita is to invert blk_nbr_3 and multiply by blk_nbr err. 
Relaxation methods and Gaussian elimination with back substitution are other means 


to the same end. 
Singular Jacobian 


If blk_nbr_j is singular, bl1k_nbr_dita can not be calculated and in the present 
incamation of WAVESIM, the simulation fails. Future versions should include an 


algorithm for attempting to recover from the singular jacobian. 
3.4.3.2.2.2.7 Correct Variables 


Each of the system variables associated with the block are corrected by 
subtracting the appropriate rows of blk_nbr dita divided by the corresponding 


element of the scaling factor vectors (sys_pot_scale or sys_flow_scale). 
3.4.3.2.2.2.8 Variable Correction Criteria 


If the block is nonlinear (blk _nbr linear flag == 0) and the variable 
correction flag is set (sb_check_var_err == 1) then blk_nbr_ive contains the 
indexes of blk_nbr_ dita which exceed in magnitude blk_nbr_imax _varcor if 
blk_nbr_alpha <1 of blk_nbr_max varcor if blk_nbr alpha 2 1. If 
blk_nbr_ive is not empty, then one of the variable corrections was too large and 
another iteration is necessary. In any case, the continuation parameter loop is 


repeated. 
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3.4.3.2.2.3 Truncation Error Control 


Once a block has been solved, a truncation error check must be performed on 
each of the associated system variables. The truncation error is assumed negligible if 
the waveform content of the last sb_nbr_we coefficients of each waveform is less than 
the limit specifed by sb_max we. If all the system variables have negligible truncation 
error, block nbr has been solved and the next block is processed. If the truncation error 
of any of the variables 1s too large, converge failure is Set to 1 to indicate the block 


has not been solved. 
3.4.3.2.3 Time Step Control: Successful Convergence 


If all the blocks successfully obtained a solution then the variable 
converge failure will equal 0. The task now is to save all of the variables in the 


history arrays, update tt 0 and tt1, update n, and update the device states. 


Update History Variables 


The history variables are extended by one column. The variable his_col 1s 


incremented and is the column index for all but the state arrays. In particular: 


his _t(1,his_col) = tt0 
his t(2,his_col) = ttl 


his N(1,his col) =N 


his nd_nbr(1:N,his_col) = var_nd_nbr 


his _ fv_name(1:N,his_col) = var_fv_name 
his _s_name(:,his_col+l) = dev_sl_name 
dev_sOQ_ name = dev_sl_name 


where (1:N,his_col) refers to the first N rows of column his_col and (:,his_col+1) 


refers to all the rows of column his_col + 1. 


Update Time Interval and Number of Coefficients 


The time interval is updated by: 
660°] ttl 


If tt0 2t1 then the simulation has successfully completed and the time loop is 


exited. Otherwise must update tt1 as well. Initially: 


ttl = ttl + ddt 


mr Dp 








Next, check if a break point (element of sb_bp) exists between tt0O and ttl. If 


such a break point exists, set tt1 equal to the earliest break point after tt0. 


Since reducing N is normally beneficial, if tt1 - tt0 > sb _dt_optimum and 
N > sb N min the algorithm assumes the waveforms are well behaved and 


decrementing N (as long aS N > sb_n_min) is appropriate. 


Since the series converged for the previous time increment, setting ddt equal to the 


minimum of 2xddt and sb_dt_max allows the system to increase the next time interval. 
Plot Intermediate Results 


Before proceeding to solve the system over the updated time interval, WAVESIM 


creates a plot of the system variables over the previous time interval. 
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3.4.3.2.4 Time Step Control: Unsuccessful Convergence 
Fatal Error 


If one of the blocks failed to converge, tt1 - tt2 < sb_dt_min, and 
N 2 sb_n_max then the simulation has failed completely and can not proceed further. In 


this case, the simulation comes to a halt prematurely. 


Recommended Recalculation Time 


If one of the blocks failed to converge and trec < tti,then tt1 = txec and the 


time loop is repeated. 


Time Increment / Number of Coefficient Control 


If one of the blocks failed to converge and  trec 2 ttl, 
ttl - tt0O < sb _dt_ optimum and N < sb_n max, N is incremented in an attempt to 
improve convergence. To improve convergence if ttl - ttO > sb dt optimum or 
N > sb _n max, the time interval is halved by setting ttl = (ttl + tt0)/2.0. 
Halving ddt is also prudent as long as ddt = sb_dt_min. Once ddt and n have been 


updated, the time loop is continued. 
3.4.3.3 Simulation Wrap-up 


Once the simulation has completed, the variables stored in the history arrays are 
plotted and saved as the user desires. If the operator desires, the device state variables 
can be used as the initial conditions for a following simulation or saved in file for future 


simulations. 


On 





3.5 Device Modelling Techniques 


The previous sections described the method WAVESIM_ uses to generate a 
mathematical system of equations and variables for interconnecting a number of different 
devices. Up to now a device has been treated as a black box characterized by its definition, 
initialization, variables which must be provided to it as resources and variables which are 


generated by it as products. As a review, here are properties of the black box: 
Definition (device. def) 


Name of Device Type 


Number of Parameters 
Names of Parameters 


Default Values of Parameters 


Number of States 
Names of States 
Default Values of State Initial Conditions 


Number of terminals 

Terminal Definitions 
Terminal Name 
Terminal Type (normal or information) 
Flow Variable Type (import or export) 
Potential Variable Type (import or export) 
Terminal KCL Group Number 


Device Structural Jacobian 


Initialization (WA VESIM input file) 


Name of Device 

Name of defining Device Type 
Parameter Values 

State Initial Conditions 


Assignment of terminals to nodes 


SO5i= 





Resources (Arguments of MATLAB device.mn file) 


Waveform type 
Import Variable Waveforms 
Parameter Values 
Value of states at beginning of time interval 
Time Structure 
Beginning time of Interval 
Ending time of interval 
Minimum time interval of interest 


Continuation Parameter 


Products (Products of MATLAB device.n file) 


Export Variable Waveforms 
Device Jacobian Matrix 
Value of states at end of time interval 
Recommended Time Structure 
Recommended Recalculation Time this interval 


Recommended ending time of next interval. 


While these specifications are the hard requirements for developing a new device type, 
they are not very constraining and it 1s possible to generate very inefficient and unworkable 


devices. The following sections are meant as guidance for developing new device types. 
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3.5.1 Import and Export Variable definitions 


One of the first tasks in designing a new object is determining which variables should 
be import variables and which should be export variables. The requirement is simply that 
the total number of export variables associated with normal terminals must equal the total 
number of import variables associated with normal terminals. To minimize the number of 
system equations however, one should usually try to define flow variables as export 


variables and potential variables as import variables. 

The constitutive equations defining a device may preclude defining all the flow 
variables as export variables. An ideal voltage source of magnitude V, for example, has 
the following constitutive equations: 

VY, = V; a5 V. 
=-f, 
Clearly, this set of equations can not be reorganized to specify both currents (flows) 


explicitly. In this case potential V, and flow J, are export variables and potential V, and 


flow J, are import variables. 
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3.5.2 Interface Variable Units 


When developing devices, a consistent convention for interface variable units is 
required. Flows are usually referenced such that positive flow into a terminal with a 
positive potential refers to power dissipated by the device. This definition is clear if the 
flow corresponds to currents or forces, but is less clear for torques. For rotating shafts 
where torques are the flow variable and rotational speed the potential, the positive 
direction for speed is in the normal operating direction while the direction for torque is 
determined by the power dissipation rule. A motor connected to a propeller would 
normally have associated a positive rotational speed and a negative torque. The propeller 
would have a positive rotational speed and a positive torque associated with its interaction 
with the motor along with a positive forward speed and negative force associated with its 
interaction with the ship dynamics. The ship dynamics model would have an associated 


positive force and positive forward speed. 


Many power system simulations go through great effort to normalize all variables by 
dividing by device base quantities to improve numerical accuracy. The models are all 
expressed in a Per Unit (PU) basis where the base quantities are machine ratings. The 
problems occur when several devices with different base quantities are combined. The 
system variables must all be scaled appropriately to ensure the elements of the system 
equations are all in the same units. Keeping the bases consistent requires much effort and 


is very prone to error. 


In WAVESIM, physical quantities using the metric system (SI) are recommended for 
all interface variables. Strict use of the metric system ensures the proper quantities are 
added and subtracted on the systems level. Individual devices may then scale the interface 
variables by their own base quantites for internal calculations. Likewise, each node of the 
system can have a scaling factor assigned to it for both flow and potential variables. In this 
manner, the beneficial aspects of the per unit system can be retained with little confusion 


as to ensuring consistent base quantities. 
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Metric System (SI) 


Length meters 

Time seconds 

Mass kilograms 
Voltage volts 

Current amperes 

Force newtons 
Angle radians 

Speed meters/second 
Rotational Speed radians/second 
Torque newton-meters 
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3.5.3 Potential References 


The node potentials are all referenced to an arbitrary value called 0. The reference 
frame for this level is a property of the device definition, but must be consistent with the 
reference frame for other device definitions to which the device may be connected. 


Following are suggested reference points: 
Electrical Voltage Volts above Ground Potential 
Mechanical Angle Radians relative to the positive vertical 


Mechanical Rotational Radians per Second relative to stationary 


Speed 

Mechanical System Dependent 

Displacement 

Mechanical Speed meters per second relative to stationary 


If mechanical rotational speeds or mechanical speeds are specified, but the actual angle is 
required within the device calculations, the speed can be integrated. If more than one 
device requires the integration of the speed, then the system modeller must ensure the state 


initial conditions corresponding to the angle or displacement is consistent for all devices. 


If an absolute reference cannot be established for a device, two terminals can be 
defined such that al] constitutive relations depend only on the difference between the two 
terminal potentials. This relative definition of potentials is commonly used for modelling 
circuit elements. An ideal transformer for example, is a four terminal device with the 
following constitutive equations: 

Vip= n(V;, Gr aaa 
he Ce hee 
lp = Lies 
lr, = = Tay 
Note that V,, is defined relative to V,,, and is a function of (V;, - V;,,). None of the export 


flow variables is a function of the absolute value of any of the potentials. 
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3.5.4 Discontinuity Control 


One of the difficulties with using vectors of orthogonal series coefficients to 
represent waveforms is the poor truncation error performance when approximating 
discontinuous variables or variables having discontinuous derivatives. Dhese 
discontinuities are usually a function of either time or the zero crossing of one of the 
variables. In any case, the time of the discontinuity is often easily determined by the 
device object. If the frequency of the discontinuities is low enough, it would be prudent 
for the device to specify the earliest discontinuity of the interval as a recommended 


recalculation time. 


If the discontinuity is a function of a waveform zero crossing, special care must be 
taken to ensure the device does not continously estimate the zero crossing to be within a 
small increment of ttO or tt1 and force the time loop to iterate tt1 around the 
discontinuity. One way around this problem is for the device to move or remove any 
discontinuities within sb_dt_ave of either tt0 or tt1 in any of its export variables. If 
sb_dt_ave is small enough, then moving the discontinuity should not affect the accuracy 
of the simulation very much yet still prevent the system time loop from hunting for the 
discontinuity by varying tt1. 

If many discontinuities occur in an export variable more frequently than sb_dt_ave, 
then the export variable should be smoothed. The smoothing operation calculates the local 
average of a waveform over the interval [t-sb_dt_ave,t+sb_dt ave]. In this manner, the 
higher order terms of the export variable are attenuated and the waveform is more likely to 


pass the truncation error test. 
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3.5.5 Consistent Initial Conditions 


Most simulation environments require the user to specify the initial values for all the 
States at time tO. In this regard WAVESIM is no different. Unfortunately, determining a 
consistent set of initial conditions which meet some definition of normal operating 
conditions 1s not an easy task for either a system modeler or a computer program. First of 
all, the concept of a normal operating condition, is not always easy to describe 
mathematically. Furthermore, even if a definition for normal operating condition, can be 


made, there is often much difficulty in determing that condition. 


An ideal solution would be for each device to calculate its own initial conditions 
during the first time increment. If a device 1s capable of determining an initial condition 
based only on its parameters and the values of its import variables, then the following 


technique can be used: 
1. Define a state called 1c, always initialize it to 0. 
2. Define Sufficient Parameters to determine the normal operating condition. 


3. Within the constitutive equations, have a check for the initial value of Ic 
equalling zero. If Ic = 0 at the beginning of the interval then use the 
equations for the normal operationg condition to determine the initial values of 
the other states. Otherwise use the initial values of the other states as passed to 


the device. In any case, the final value for the state 1c should be set to 1. 


This method for determining the initial conditions is well suited for determining the 
initial conditions of the states of rotating machines. Essentially, a load flow is conducted 


in the first tinne increment to determine the initial state values. 
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3.5.6 Waveform type conversion 


Performing the calculations for the constitutive equations for certain devices may be 
easier to accomplish in one waveform over other waveforms. Converting the import 
variables to a fixed waveform type is permissible and at times desirable. As long as the 
export variables are converted back to the proper type and the jacobians reflect the 


waveform conversions, all should work out well. 


If the export variables depend on higher order terms of intermediate calculations, 
converting the import variables to waveforms of a length longer than w and performing all 
of the intermediate calculations using this longer length before truncating back to N when 


generating the export variables may be desirable in avoiding excessive truncation errors. 
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Chapter 4 WAVESIM 
4.1 Basic Description 


WAVESIM, a simulation program written in the C programming language, 
demonstrates the algomthms discussed in detail in Chapter 3 for simulating systems of 
nonlinear lumped parameter models representing the electro-mechanical components 
comprising an Integrated Electric Drive system. The general characteristics of WAVESIM 
are: 

1. | System and Simulation Parameters specified in a text Input File. 


Zz Device Definitions are in text file device. def. 


3. WAVESIM Performs following 4 tasks: 
A. Reads in Device Definitions and initializes simulation. 
B. Reads Input File and determines devices and nodes of system. 
C. Builds and reduces system into a sequence of blocks. 
D. Writes a MATLAB Script file for conducting the simulation. 


4. The actual Simulation is conducted in MATLAB. 


5. Supported Waveform types are: 
Data Series. 

Fourier Series. 

Legendre Series. 


Polynomial Expansions. 


moO wD > 


Chebyshev Series. 
6. | Waveform operators are MATLAB functions defined in M-files. 


7. Device Constitutive Equations are detailed in MATLAB functions defined in 
M-files. 


8. The present Incarnation of WAVESIM has these limitations: 
A. Subsystems have not been implemented. 
B. System and Device Structural Jacobians must be time independent. 
C. Newton-Raphson is the only equation solving method used. 


Relaxation Techniques have not been implemented. 
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MATLAB was chosen as the environment for conducting the simulation for the 


following reasons: 


Ne 


MATLAB is ideally suited for treating vectors and matrices as abstract data 
types. 

MATLAB has built in plotting routines. 

The ability to create MATLAB M-files which when invoked, execute a long 


series of commands called a script. M-files can also be used to create new 
MATLAB functions. 


MATLAB has many built in functions for analysing matrix properties. 


Since WAVESIM is an algorithm demonstration program, speed is not of 
primary concem. Interest in determining if the algorithms work is of higher 


interest than optimizing for speed. 
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4.2 Running WAVESIM 


Under either the UNIX operating system or IBM DOS, WAVESIM is executed by 


entering at the commmand prompt: 
athena% wavesim file.in 
where £ile.in is an optional entry for the file name of the input file. WAVESIM 


will attempt to read in the device.def file and if successful, will display the following 


header: 


WAVESIM 
Revision 2.0 <> Ro pap ila eo Ae 


(Cj Copyright 1990,1991 by Norbert H. Deerry 


If WAVESIM encountered errors when reading device.def, an error message is 


printed and the program terminates. 


If file.in was not specified on the command line, the user is prompted for a file 


Name. 


Enter WAVESIM INPUT file name 


If instead of a file name q is entered, WAVESIM terminates execution. A directory 
listing can be obtained by entering a ? followed optionally by a file specification (operating 
system dependent). 


Under normal execution of WAVESIM, there is no further interaction with the user. 
WAVESIM automatically creates an output file having the same base filename as 


file.in but having .m as an extension (i.e. file.in becomes file.m). 


NOTE: Do not create input files with .m extensions as these files will be overwritten 
by WAVESIM. Also avoid using file names which are valid MATLAB functions. 


WAVESIM provides extensive support for providing the user with feedback through 
the use of the DEBUG command. Most of the major routines in WAVESIM have a debug 
option for displaying the results of calculations internal to WAVESIM. 
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If errors are found reading either device. def or the input file, WAVESIM displays 
an error message which includes the file name and the line number within the file. 
WAVESIM attempts to continue reading an input file even if errors are detected but will 


only create an output file if no errors are encountered. 
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4.3 Input File Specification 


The Input File describes the system topology, defines the device parameters, and 


specifies simulation paramaters. The basic charactenstics of the file are: 


1. 


2 


mee 


ASC II text files. 
Lines beginning with %, # or ! are ignored. Empty lines are ignored as well. 


Data lines can be continued on the following line if the last characters in the line 


ANG 7. OF \" 


Commands all begin with a key-word. Key-words are case insensitive and 
usually can be truncated to three letters unless a conflict with another key-word 


exists. 

Commands and their arguments may be separated by either spaces or tabs. 

The contents of other files can be incorporated by using the INCLUDE command. 
Single Line Commands have data arguments entered on only one line. 


Multiple Line Commands consist of groups of subordinate commands. The group 


must end with a line beginning with the key-word END. 


Here is a summary of the Commands available : 


DEBUG Print Debug Information 

DEFAULT Default System Parameter Initialization 
DEVICE Device Specification 

INCLUDE Include another file 

NODE Node Parameter Specification 

TIME Time Increment Control 
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% 

eeceooc , Ln 

% 

debug 
busild» system identity 
build system_blocks 
fanned block 
END 

% 


device VDC SOURCE Vs 
TERMINAL 1 1 
TERMINAL 2 0 
PARAMETER VS 1.0 
END 

% 

device RESISTOR Ri 
TERMINAL 11 
TERMINAL 2 2 
PARAMETER Re WOR) 
END 

% 

device RESISTOR R2 
TERMINAL 1 2 
TERMINAL 2 3 
PARAMETER Re. 160 
END 

% 

device INDUCTOR Li 
TERMINAL 1 2 
TERMINAL 2 0 
PARAMETER io 6 
END 

% 

device INDUCTOR L2 
TERMINAL 1 3 
TERMINAL 2 0 
PARAMTER 1 eee even 6, 
END 

% 

device CAPACITOR Cl 
TERMINAL 1 2 
TERMINAL 2 0 
PARAMETER Co 0 
END 

% 

device CAPACITOR C2 
TERMINAL 1 3 
TERMINAL 2 0 
PARAMETER Ce. BiG 
END 


Example Input File 


=< 
" 
| 


}+— 


© 


L=ld 


L=0 | 


a 


Jad 


l=¢J 
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5 
node 


% 


a 
scale 
scale 
error 
error 
end 


default 


time 


Gmin 
Rmin 


check 


SCLlLrEOrE 
SCrErEOr 
LS as a @ 2 
CrcOor 
error 
error 
Sit a Gi @ B a 
error 


scale 
scale 


Example Input File (continued) 


Bot ent ial. 1.0 
flow a>.'0 
kel 5e-3 
pot 5e-s 


0 
0 


both 


kel 
pot 
node 
flow 
kel 
Doe 
node 
flow 


eqn 
eqn 
Var 
var 
mult 
mult 
WL 
mult 


potential 1. 
flow i 


max count. 10 
max int count 6 


alpha 
alpha 
alpha 


my lan © 
slp gh emma 9 Wolk ot 
ric min 


eZo 
ry Os. 


diverge start 3 
diverge max cnt 2 
diverge error mult 10.0 


waveform 
waveform 


wtype 
nbr 
nbr 
nbr 
nbr 


END 


start 


finish, Z 


END 


min) 0 
max 5 
Ope. 0 
eM iete: als, 
ave QO 

0 

0 


content max .005 
Content nbE Z 


3 


coef 7 
coef 
coef 
data 


min 6 
max 14 
20 


2025 


2250 


a 
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4.3.1 DEBUG 


If DEBUG 1s specified without any arguments, the command is interpreted as a 
multi-line command. Each of the following lines should contain the name of one of the 
subroutines listed below. If the key-word OFF follows the subroutine name, the debug flag 
for that subroutine is turned off. Otherwise, the debug flag for the specified routine 1s 


turned on. The last line of the group should begin with the key-word END. 


If DEBUG is specified with arguments, the command is interpreted as a single-line 
command and the arguments should consist of one of the subroutines listed below and 
optionally, the key-word OFF. A single line command does not have an END keyword 


associated with it. 


Here is a list of subroutines for which debug flags have been defined (Note: The 


subroutine names are case sensitive) 


init devices 

read device def 
read file 

read file device 
read file default 
read file node 

read file time 

read file debug 

build system 

build system_identify 
build system _structural_jacobian 
build system _blocks 
find block 

print system identify 


write file 
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4.3.2 DEFAULT 


If DEFAULT is specified without any arguments, the command is interpreted as a 
multi-line command. Each of the following lines should contain one of the subordinate 


commands listed below. The last line of the group should begin with the key-word END. 


If DEFAULT is specified with arguments, the command is interpreted as a single-line 
command and the arguments should consist of one of the subordinate commands listed 


below. A single line command does not have an END keyword associated with it. 


Here is a summary of the DEFAULT subordinate commands: 


ALPHA Continuation Parameter Control 
CHECK Error Checking Flags 

DIVERGE Divergence Test Control 

ERROR Default Error Levels 

GMIN Default Node Leakage Conductance 
MAX Maximum Iteration Counts 

NBR Number of Coefficients Control 
RMIN Default Node Series Resistances 
SCALE Default Variable Scaling Factors 


WAVEFORM CONTENT 


WTYPE 


Waveform Content Limits 


Waveform Type 
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4.3.2.1 DEFAULT: ALPHA 


The ALPHA subordinate command specifies the parameters needed to control the 


continuation parameter for nonlinear blocks. 


Command Description MATLAB Variable 
ALPHA INIT Value Continuation Parameter sb_alpha_init 
Initial Value 
ALPHA INC INIT Value Continuation Parameter sb_dalpha_init 


Initial Increment 


ALPHA INC MIN Value Minimum Continuation sb_dalpha_min 
Parameter Increment 


For a nonlinear block, the continuation parameter is initialized to the ALPHA INIT 
value. The initial increment for the continuation parameter is specified by ALPHA INC 
INIT. If the block fails to converge, the continuation parameter 1s progressively 
decremented until the block converges or if convergence fails due to the difference 
between the last value of the continuation parameter that converged and the present value 
of the continuation parameter being less than ALPHA INC MIN. If the block converges, 
the continuation parameter 1s incremented by ALPHA INC INIT until it equals I. 


4.3.2.2 DEFAULT: CHECK 


The CHECK subordinate command determines for nonlinear blocks, whether the 
equation error, the variable correction magnitude, or both should be used for the 


convergence criteria. 


Command Description MATLAB Variable 

CHECK EQN Check only Equation Errors sb check_eqn_err = 1 
sb check var _err = 

CHECK VAR Check only Variable sb check eqn err = 0 

Corrections sb check var err = 1 

CHECK BOTH Check both Equation Errors sb_check_eqn_ err = 1 

and Variable Corrections sb check var err = 1 
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4.3.2.3 DEFAULT: DIVERGE 


The DIVERGE subordinate command specifies when and how to check a nonlinear 
block for divergence. After DIVERGE START iterations, if the largest relative error 
increases for DIVERGE MAX CNT iterations and the relative error is at least DIVERGE 
ERROR MULT then the block is assumed to be diverging and the converge_failure flag 


iS set. 


Command Description MATLAB Variable 


DIVERGE START Value Number of iterations to wait sb_div_start_cnt 
before testing for divergence 


DIVERGE MAX CNT Value Number of iterations to sb div_max_cnt 
allow relative error. to 
increase befor concluding 
divergence 


DIVERGE ERR MULT Value Value of relative error sb i div_err 
below which to ignore 
divergence iteration count 


Seve 








4.3.2.4 DEFAULT: ERROR 


The ERROR subordinate command determines for nonlinear blocks, the default 
maximum equation errors and variable corrections which are permissible. These default 
values can be overridden for a specific node with the NODE command. When the 
conituation parameter equals 1, ERROR EQN KCL is the maximum error for the node 
KCL equtions and ERROR EQN POT is the maximum error for the potential difference 
equations. Likewise when the continuation parameter equals 1, ERROR VAR NODE is 
the maximum correction to a node potential and ERROR VAR FLOW is the maximum 
correction to an import flow variable. The ERROR MULT subordinate commands are 


multipliers to the above limits for continuation parameters less than I. 


Command Description MATLAB Variable 

ERROR EQN KCL Value Default KCL Equation sys_kel_err’ 
Maximum Error 

ERROR EQN POT Value Default Potential Difference sys _var_err’ 
Maximum Error 

ERROR VAR NODE Value Default Maximum sys_nd_err’ 
correction to Node 
Potentials 

ERROR VAR FLOW Value Default Maximum sys_fv_err’ 
correction to Import Flow 
Variables 

ERROR MULT KCL Value Multiplier to ERROR EQN sb i kel err 


KCL when continuation 
parameter < 1 


ERROR MULT POT Value Multiplier to ERROR EQN sb_i_ pot err 
POT when _— continuation 
parameter < 1 


ERROR MULT NODE Value Multiplier to ERROR VAR sb_i nd err 
NODE when_ continuation 
parameter < ]| 


ERROR MULT FLOW Value Multiplier to ERROR VAR sb i fv_err 
FLOW when continuation 
parameter < | 


Note 1: sys_xxx_err are actually arrays containing for each equation or variable, either 
the default value specified here or the overriding value specified in the NODE 
command. 


-115- 





4.3.2.5 DEFAULT: GMIN 


The GMIN subordinate command defines the default value for G,,,. G,,;, 18 used to 
modify the KCL equations to help prevent singular systems. G,,,, Should normally be set 


to 9 unless a singularity problem exists. The value for G_,, can be overridden for a 


R 


particular node through the NODE command. 


Command Description MATLAB Variable 


GMIN Value Leakage Conductance to 0 sys _Gmin’ 
Potential 


Note 1: sys Gmin 1s actually an array containing the value for G,,, for each node: either 
the default value specified here or the overriding value specified in the NODE 
command. 


4.3.2.6 DEFAULT: MAX 


The MAX subordinate command determines the maximum number of 
Newton-Raphson iterations for a nonlinear block before the continuation parameter 1s 
decremented. MAX COUNT specifies the maximum number of iterations when the 
continuation parameter equals 1 while MAX INT COUNT specifies the maximum 


number of iterations when the continuation parameter 1s less than 1. 


Command Description MATLAB Variable 
MAX COUNT Value Maximum number — of sb_maxcnt 
Iterations when the 
continuation parameter 
equals 1 
MAX INT COUNT Value Maximum number of sb_i _maxcnt 
Iterations when the 
continuation parameter 1s 
less than 1 
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4.3.2.7 DEFAULT: NBR 


The NBR subordinate command controls the number of coefficients the waveforms 
will have. NBR COEF specifies the initial number of coefficients to use. NBR COEF 
MIN is the minimum number of coefficients to use while NBR COEF MAX 1s the 
maximum number of coefficients). NBR DATA is the number of data points per 


waveform used when generating plots. 


Command Description MATLAB Variable 

NBR COEF Value Initial number of n 
coefficinets 

NBR COEF MIN Value Minimum number of sb_n min 
coefficients 

NBR COEF MAX Value Maximum number of sb_n max 
coefficients 

NBR DATA Value Number of points per sb_n data 


waveform to use in plots. “ 


4.3.2.8 DEFAULT: RMIN 


The RMIN subordinate command defines the default value for R,;,. Rmmin 1S used to 
modify the Potential Difference equations to help prevent singular systems. R,,,, should 
normally be set to 0 unless a singularity problem exists. The value for R,,;, can be 
overridden for a particular node through the NODE command. 


Command Description MATLAB Variable 
RMIN Value Series Resistance for Export sys Rmin’ 
Potentials 


Note 1: sys _Rmin is actually an array containing the value for R,,;, for each node: either 
the default value specified here or the overriding value specified in the NODE 
command. 
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4.3.2.9 DEFAULT: SCALE 


The SCALE subordinate command specifies the default scaling parameters for the 
potential and flow variables. The default scaling parameters can be overridden for a 


particular node through the NODE command. 


Command Description MATLAB Variable 
SCALE POTENTIAL Value Default scaling factor for sys pot scale’ 
Potentials 
SCALE FLOW Value Default scaling factor for sys flow scale’ 


Flow Variables 


Note 1: sys pot scale and sys flow scale are actually arrays containing the scaling 
factors for each node: either the default values specified here or the overriding 
values specified in the NODE command. 


4.3.2.10 DEFAULT: WAVEFORM CONTENT 


The WAVEFORM CONTENT subordinate command controls the maximum 
allowable truncation error by specifying the maximum waveform content WAVE CONT 
MAX for the last WAVE CONT NBR coefficients of a waveform. 


Command Description MATLAB Variable 
WAVE CONT MAX Value Maximum Waveform sb_max hh 
Content 
WAVE CONT NBR Value Number of Coefficients to sb nbr_hh 


apply maximum to. 


4.3.2.1] DEFAULT: WTYPE 


The WTYPE subordinate command specifies the waveform type to use in the 


simulation 
Command Description MATLAB Variable 
WTYPE Value Waveform Type Indicator wtype 


Data Series 

Fourier Series 

Legencre Series 
Polynomials 

MATLAB Polynomials 
Chebyshev Series 


Nn B&B WN = 
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4.3.3 DEVICE 


DEVICE is always a multi-line command. The command must be entered in the 


following format: 


DEVICE Device Type Name 


where: 
Device _Type Device Type Name from device.def file. 


Name Name of this particular device. 


The subordinate commands for the DEVICE command are: 


TERMINAL Assign Terminals to Nodes (mandatory). 
PARAMETER Assign Parameter Values (optional). 
STATE Assign State Initial Conditions (optional). 


All of the terminals as defined in the device.def must be assigned to a node. If the 
parameters or states are not assigned values, the default values specified in device .def 


are used. 


The last line of the command group must begin with the key-word END 


4.3.3.1 DEVICE: TERMINAL 


The TERMINAL subordinate command assigns a terminal to a node and must be 


entered in the following format: 


TERMINAL Terminal Name Node Nbr 


where: 
Terminal Name Terminal Name from device.def file. 


Node Nbr Serial Number of Node this terminal is attached to. 


All of the terminals as defined in device .def must be attached to a node of the 


system. 
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4.3.3.2 DEVICE: PARAMETER 


The PARAMETER subordinate command assigns a value to a parameter of the 
device. If the parameter is a single value as defined in device. def then the parameter 


command must be of the following format: 


PARAMETER Paramet er Name Value 


where: 
Parameter Name Parameter Name from device. def file. 
Value Parameter Value. 
If the parameter 1s a matrix as defined in device.def then the parameter 
command must be of the following format: 


PARAMETER Parameter Name MATRIX 
ma trix val ues 


END 


where: 


Parameter Name Parameter Name from device. def file. 


matrix values Parameter matrix. The number of rows and columns 
of the matrix must be the same as specified in the 
device .def file. Rows are entered one line at a 
time with columns separated by spaces. 
If a parameter as defined in device. def is not assigned a value, then the default 


values specified in device. def is used. 
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4.3.3.3 DEVICE: STATE 


The STATE subordinate command assigns an initial value to a state of the device 


and must be entered in the following format: 


STATE State Name Value 


where: 
State Name State Name from device.def file. 


Value Initial value of state. 


If a state as defined in device.def is not assigned a value, then the default 


values specified in device. def is used. 
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4.3.4 INCLUDE 


INCLUDE is always a single-line command. The command must be entered in the 


following format: 


INCLUDE File Name 
where: 
File Name Name of the file to include. 


The contents of the included file are inserted at the location of the INCLUDE 


command. 
4.3.5 NODE 


NODE is always a multi-line command. The command must be entered in the 


following format: 


NODE Node Nbr 


where: 


Node Nbr Serial Number of the node. 


The subordinate commands for the NODE command are: 


ERROR Node Error Levels 

GMIN Specify node G,,;, Value 
NAME Assign a name to the node 
RMIN Specify node R,,;, Value 
SCALE Specify node scaling factors 


The last line of the command group must begin with the key-word END 
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4.3.5.1 NODE: ERROR 


The ERROR subordinate command determines for nonlinear blocks, the maximum 
equation errors and variable corrections which are permissible. These values override the 
default values. When the conituation parameter equals 1, ERROR EQN KCL is the 
maximum error for the node KCL eqution and ERROR EQN POT is the maximum error 
for the potential difference equations. Likewise when the continuation parameter equals 
1, ERROR VAR NODE is the maximum correction the node potential and ERROR VAR 
FLOW is the maximum correction to an import flow variable. The ERROR MULT 
subordinate commands of the DEFAULT command are multipliers to the above limits for 


continuation parameters less than 1. 


Command Description MATLAB Variable 
ERROR EQN KCL Value Maximum KCL _ Equation sys_kcl_err’ 
Error 
ERROR EQN POT Value Maximum Potential sys _var_err’ 


Difference Error 


ERROR VAR NODE Value Maximum correction to sys_nd_err’ 
Node Potential 


ERROR VAR FLOW Value Maximum correction to sys fv_err’ 
Import Flow _ Variables 
attached to this node 


Note 1: sys_xxx_err are actually arrays containing for each equation or variable, either 
the default value or the overriding value specified here. 


4.3.5.2 NODE: GMIN 


The GMIN subordinate command defines the value for G,,,,.._ G,,;, 18 used to modify 
the KCL equation to help prevent singular systems. G,,;,, should normally be set to 0 


unless a singularity problem exists. The value for G,,,, overrides the default value. 


Command Description MATLAB Variable 


GMIN Value Leakage Conductance to 0 sys Gmin’ 
Potential 


Note 1: sys _Gmin 1s actually an array containing the value for G,,,, for each node: either 
the default value or the overriding value specified here. 
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4.3.5.3 NODE: NAME 


The NAME subordinate command 
Command Description MATLAB Variable 
NAME Node Name Name of the Node sys node name’ 


The node name is only used to associate the node serial number to a more 
understandable label. The node name is optional and does not affect computation in any 


way. 
Note 1: sys_node_name Is actually an array containing the names of all the nodes. 


4.3.5.4 NODE: RMIN 


The RMIN subordinate command defines the node value for R,,,. R,;, 1S used to 
modify the Potential Difference equations to help prevent singular systems. R,,,,, Should 
normally be set to 0 unless a singularity problem exists. The value for R,;, overrides the 
default value. 


Command Description MATLAB Variable 
RMIN Value Series Resistance for Export sys Rmin’ 
Potentials 


Note 1: sys Rmin is actually an array containing the value for R,,,, for each node: either 
the default value or the overriding value specified here. 


4.3.5.5 NODE: SCALE 


The SCALE subordinate command specifies the node scaling parameters for the 


potential and flow variables. The scaling parameters override the default values. 


Command Description MATLAB Variable 
SCALE POTENTIAL Value Scaling factor for Node sys pot scale’ 
Potential 
SCALE FLOW Value Node scaling factor for Flow sys flow scale’ 
Variables 


Note I: sys pot scale and sys flow scale are actually arrays containing the scaling 
factors for each node: either the default values or the overriding values specified 
here. 
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4.3.6 TIME 


If TIME is specified without any arguments, the command is interpreted as a 
multi-line command. Each of the following lines should contain one of the subordinate 


commands listed below. The last line for the section should begin with the key-word END. 


If TIME is specified with arguments, the arguments should consist of one of the 
subordinate commands listed below. A single line command does not have an END 


keyword associated with it. 


The subordinate commands for the TIME command are: 


BREAK Insert Break Point 

DT Time Increment Control 
FINISH Ending Time of Simulation 
START Starting Time of Simulation 


4.3.6.1 TIME: BREAK 


The BREAK subordinate command inserts a simulation break point which forces a 
waveform boundary to occur at the designated time. Bracketing intervals in which a 
discontinuity will occur with breakpoints can reduce the computational effort required by 
WAVESIM. 


Command Description MATLAB Variable 
BREAK Time Break Point time sb_ bp’ 
sb bp _nbr 


Note 1: sb_bp is actually an array of break points in chronological order. sb _bp_nbr is 
the number of break points. 
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4.3.6.2 TIME: DT 


The DT subordinate command controls the wavefonn interval. 


Command Description MATLAB Variable 

DT MIN Value Minimum Waveform sb _dt_min 
Increment 

DT MAX Value Maximum Waveform sb _dt_max 
Increment 

DT OPTIMUM Value Optimum Waveform sb _dt_ optimum 
Increment 

DT INITIAL Value Initial Waveform Increment sb_dt_ init 

DT AVE Value Minimum Time Interval of sb_dt_ave 


Interest (Averaging Interval) 


If the time interval is less than DT OPTIMUM, the number of coefficients is less 
then the maximum and a block does not converge, the number of coefficients 1s increased 
for the next iteration. Otherwise, if the block does not converge the time interval is 


reduced. 


DT AVE is the minimum time interval of interest and is used by devices to smooth 


their export waveforms or to move discontinuity boundaries. 


4.3.6.3 TIME: FINISH 


The FINISH subordinate command specifies the ending time of the simulation. 


Command Description MATLAB Variable 


TIME FINISH Value Ending Time of Simulation t1 


4.3.6.4 TIME: START 


The START subordinate command specifies the beginning time of the simulation. 


Command Description MATLAB Variable 


TIME START Value Starting Time of Simulation t0 
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4.4 Device Definition File Specification 


The device definition file device.def contains the definitions of the device types 
which can be specified in a WAVESIM Input File. The basic characteristics of the file are: 


1 
Z. 


oe 


ASC II text files 
Lines beginning with %, # or ! are ignored. Empty lines are ignored as well. 


Data lines can be continued on the following line if the last characters in the line 


Ave “42 OF \. 


Commands all begin with a key-word. Key-words are case insensitive and 
usually can be truncated to three letters unless a conflict with another key-word 


exists: 

Commands and their arguments may be separated by either spaces or tabs. 

The contents of other files can be incorporated by using the INCLUDE command. 
Single Line Commands have data arguments entered on only one line. 


Multiple Line Commands consist of groups of subordinate commands. The group 


must end with a line beginning with the key-word END. 
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% 
: device.def 


Example device.def File 


& debug init devices 
Se meaa device deft 


device RESISTOR 


Terminal 1 
Terminal 1 
Terminal 2 
Terminal 2 


Parameter 
EPunce on 
Seructural 
DD 
end 
end 
% 


device INDUCTOR 
Terminal 1 
Terminal 1 
Terminal 2 
Terminal 2 
Parameter 
Function 
Structural 
LL 
LL 
end 
end 
% 


EOt ViT import 
Flow it Export 
Pot V2 Import 
Flow 12 Export 
Rete—is 

resistor 


Jacobian All 


POL Vi dmpore 
Flow PA ES port 
Pot V2 Import 
Flow 2 Export 
L..te=15 

inductor 


Jacobian All 


device CAPACITOR 


Terminal 1 
Terminal 1 
Terminal 2 
Terminal 2 
Parameter 
Function 

Structural 


OD 
end 
end 
% 


POt Vil Expore 
Flow i YExpOrE 
Pot V2 Import 
Flow IZ wmpore 
C le-i35 
Capacitor 


Jacobian All 


device VDC SOURCE 


Terminal 1 
Terminal 1 
Terminal 2 
Terminal 2 
Parameter 
Function 

Structural 

10 


OD 
end 
end 
% 


POL Vi Export 
Flow id EXpore 
BOL V2 Import 
Flow TZ wlmport 
VS 1.0 

VdeC=Sre 


JacObian All 


device REFERENCE 


Terminal Gnd Pot 


VO Export 


Terminal Gnd Flow aN Impex t 


Parameter 
SEructural 
0 


end 
end 


Vref 0. 
Jacobian ALL 


% include load flow definitions 
acide loadflow.def 


peanclude rotatin 
fee oe powersys. 


% include other circuit elements 


include etre .clmeder 


g machinery IED 
ef 


Models 
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4.4.1 DEBUG 


The DEBUG command is always a single-line command and results in the display of 


debug information for a specified routine during the execution of WAVESIM. 
Command Description 
DEBUG init devices Print Info on Initial System Parameters 


DEBUG read device def Print Info on what is read from device . def 


4.4.2 DEVICE 


DEVICE is always a multi-line command. The command must be entered in the 


following format: 
DEVICE Device _ Type 
where: 
Device Type Device Type Name (must be unique) 


The Device Type Name is used to correlate a given device in an input file with the 


properties of the device as specified here. The subordinate commands for DEVICE are: 


TERMINAL Specify Terminal Variable Properties 
PARAMETER Specify Parameters 

STATE Specify States 

FUNCTION Specify MATLAB function 


STRUCTURAL JACOBIAN Specify Structural Jacobian 


The last line of the command group must begin with the key-word END 
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4.4.2.1 DEVICE: TERMINAL 


The TERMINAL subordinate command defines the properties of the variables 
associated with a terminal. If the Terminal is a normal terminal, both the flow and 
potential variables need definitions. Flow variables also require a KCL group number 
KCL which corresponds to the group of terminals for which KCL can be written internally 
to the device. If the flow variable does not belong to a KCL group, its value should be 0. 
Variable are IMPORT 1s they are a resource to the device and are EXPORT if they are a 
product of the device. The total number of export variables associated with normal nodes 
must equal the total number of import variables associated with normal nodes. 


Normal Node potentials are defined by either 


TERMINAL Terminal Name POTENTIAL Variable Name EXPORT 


or 


TERMINAL Terminal Name POTENTIAL Variable Name IMPORT 


Normal Node flows are defined by either 


TERMINAL Terminal Name FLOW Variable Name EXPORT KCL 


Or 


TERMINAL Terminal Name FLOW Variable Name IMPORT KCL 


Information Node potentials are defined by either 


TERMINAL Terminal Name INFORMATION Variable Name EXPORT 


or 
TERMINAL Terminal Name INFORMATION Variable Name IMPORT 
Where 
Terminal Name One word name for Terminal 
Variable Name One word name for Variable 
KCL KCL Group Number (0 if none) 
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4.4.2.2 DEVICE: PARAMETER 


The PARAMETER subordinate command defines the parameters of the device and 
optionally, declares the default values for the parameters. Parameters can either be single 


valued or a matrix. A single valued parameter is defined by: 


PARAMETER Paramet er Name Defaul t_Val ue 


where 


Parameter Name One word name for Parameter 


Default_Value Optional Default Value for Parameter 


Matrix parameters for which which no default values are provided are defined by: 


PARAMETER Parameter Name MATRIX Nbr_Row Nbr_Col 


where 


Nbr_Row Number of rows in Matnx 
0 should never be used 
-] indicates variable dimensioned 


NprwCoi Number of columns in Matnx 
Q should never be used 
-] indicates variable dimensioned 


Matrix parameters for which which default values are provided are defined by: 


PARAMETER Parameter Name MATRIX Nbr Row Nbr Col 
DEFAULT 

Default Matrix 

END 


where 


Default Matrix Default Matrix values, Each matrix row should 
be entered one line at a time with columns 
Separated by spaces. 
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4.4.2.3 DEVICE: STATE 


The STATE subordinate command defines the states of the device and optionally, 


declares the default initial values for the states. States are defined by: 
STATE State Name Default Value 

where 
State Name One word name for State 


Default Value Optional Default Initial Value for State 
4.4.2.4 DEVICE: FUNCTION 


The FUNCTION subordinate command is mandatory and defines the MATLAB 
function which defines the device constitutive equations. The MATLAB function is 
specified by: 

FUNCTION MATLAB Function 
where 


MATLAB Function MATLAB Function name 
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4.4.2.5 DEVICE: STRUCTURAL JACOBIAN 


The STRUCTURAL JACOBIAN subordinate command defines the structural 
jacobian matrix of the device for 1 or for all of the waveform types. The structural 


jacobian for all waveform types is specified by: 


STRUCTURAL JACOBIAN ALL 
Struct ural Jacobian 


END 


where 


Structural Jacobian Structural Jacobian Matrix. The Rows 
correspond to Export Variables ordered 
according to the order of definition. Similarly, 
the Columns correspond to Import Variables 
ordered according to the order of definition. 
The elements are Structural Jacobian Codes 
detailed below. 


The structural jacobian for one particular waveform type is specified by: 


STRUCTURAL JACOBIAN Waveform Type 
Struct ural Jacobian 


END 


where 


Waveform Type Waveform Type Code the structural jacobian is 
defined for 
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Structural Jacobian Codes 


—— a 


Code {Type of Matnx 





— a ee 







Zero Matrix (all elements are always zero) 


Meat ee meaner ont 


Note 1: An AC Matrix is one for which the constant component of the export variable 
depends only on the constant component of the import variable. The other 
components of the export variable can not depend on the constant component of the 
import variable but are not restricted in any other way. 


rit 
| 
15 
i 








Waveform Type Codes 


Waveform Type 
<a 


Legendre Series 


Polynomials 


_ eT 


Matlab Polynomials 


Chebyshev Series 





eee 





4.4.3 INCLUDE 


INCLUDE is always a single-line command. The command must be entered in the 


following format: 


INCLUDE File Name 


where: 


File Name Name of the file to include. 


The contents of the included file are inserted at the location of the INCLUDE 


command. 
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4.5 Adding devices 


Adding New devices to WAVESIM requires the creation of a MATLAB M-file 
defining the device constitutive equations and the addition of an entry in the device. def 
file. 


4.5.1 MATLAB M-FILE 


Creating a MATLAB M-FILE for generating a new device requires adherence to a 
strict function argument list format. The following header indicates the format required by 
WAVESIM: 


function [e, jacob, s1,tt1])=function’ (stype,i, par,mparl’,mpar2,s0,tt, alpha) 
FUNCTION 


VERSION 2-5 of 19 April 1991 
(GC) Copyright. 1990, 1991 “by Norbert HH. Doerry 


[e , jacob, sl, ttl] = function(stype,i, par,s0,tt, alpha) 


FUNCTION creates the values and jacobian matrix for a FUNCTION 


data points 
fourier series 
legendre series 
polynomial 

MATLAB Polynomials 
chebyshev series 


SLYPpe 


it uu we al 
HUB WN HF 


[il i2 ...] where 
il, iZ, ... are column vectors of import variables 


pP- 
II 


par = [pl p2 ...] where 
parameter 1 
parameter 2 


pl 
p2 


lt il 


mparl = matrix parameter parameter Ml 


mparZ = matrix paraemter parameter M2 
s0 = [S01 S0 2 ...] where 
SO 1 = state _ 1 value at £0 
SQ 2 = state 2 value at £0 
ct = [tO tl dtave] where 
EO = initial time of the interval 
td = final time of the interval 
dtave = averaging increment 
alpha = continuation parameter 


dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP 
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% = = [el e2 ...}) where 

% el, e2, ... are column vectors of export variables 
% 

% jacob = Jacobian matrix of e with respect to i 

% 

% sl = ed edd ee || where 

% 7 Dim d eae seetee 2 ane atl 

% S12 = state 2 value at tl 

% eit 

% 

% 

% ee = [ntl ntt]}] where 

% ntl = recommended recomputation time this interval 
% ntt = recommended ending time next interval 

% 


Note 1: function is the name of the function defining the device. The MATLAB M-FILE 
should be called function. m. 


Note 2: mpar1, mpar2, etc. are only specified if the device as defined in device.def 
has matrix parameters. 
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Example MATLAB M-File 


hameceron [e@ , Jacob, Si, ttl]~= tesister(stype,i,par,s0,tt;,alipnha) 


OP dP dP dP dP OP OP OP OP OP OP dP OP dP dP dP OP DP dP dP dP DP OP oP dP OP oP dP OP oP oP OP 


2 J oP dP OP OP OP 


ct 
© 


ge Q. ct 
re 


vs wg 


% 


r= 
% 


sl 
% 
vl 
nt a 
% 
ae 


p 
NO 


OP LJ. dP H+ cP oP (D oP oP 


}- 


RESISTOR 


VERSION 1.6 of 25 February 1991 
(C) Copyright 1990,1991 by Norbert H. Doerry 


[e , jacob, sl, ttl] = resistor(stype,i,par,s0,tt,alpha) 


resistor creates the values and jacobian matrix for a resistor 


stype = 1 data points. 
= 2 fourier series 
= 3 see peer series 
= 4 polynomial 
st = [vl v2] where vl and v2 are column vectors 
Bae = an where Rais the resistance 
Ss = 
eo =P [tOrtel at.) . 
tO =wolnmatagl. €alie, Of tne. 2uterval 
GL jsiinal- time of the amterval 
at = averaging time interval 
alpha = continuation parameter 
e = [i1 12] where il and i2 are column vectors 
jacob = Jacobian matrix of e with respect to i 
sl = 
me = Oe ntt]) where 


ntl erecommencea Trecomputataon tame this interval 
ntt recommended ending time next interval 


structural jacobian 


= [tl t0); 


par (l); 


Ge 


Hcy); 
nig ae) Mes 


(vl - v2) / R; 
- il; 


[fr 22); 


= eye (n); 


feeeo-={ i[- 41 (/oR 1177, R= = 42/7 R.i1°57 R ); 
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4.5.2 device.def File 


A DEVICE entry must be made in the device. def file as described in a previous 


section. Here is an example of the entry made for the resistor: 


% 
DEVICE RESISTOR 
% 
TERMINAL 1 POTENTIAL Vl IMPORT 
TERMINAL 1 FLOW II EXPORT 1 
TERMINAL 2 POTENTIAL V2 IMPORT 
TERMINAL 2 FLOW IZ ,°EXPORT 1 
PARAMETER R le~-15 
FUNCTION resistor 
STRUCTURAL JACOBIAN ALL 
DD 
DD 
END 
END 


Note: a device can be have multiple entries in device.def to reflect different 
default state initial values and default parameter values. For example, one may desire to 


create a model of a 1000 ohm resistor: 


% 
DEVICE 1K RESISTOR 
% 
TERMINAL 1 POTENTIAL V1 IMPORT 
TERMINAL 1 FLOW Et EXPORT 1 
TERMINAL 2 POTENTIAL V2 IMPORT 
TERMINAL 2 FLOW I2 EXPORT 1 
PARAMETER R 1000 
FUNCTION resistor 
STRUCTURAL JACOBIAN ALL 
DD 
DD 
END 
END 


In this manner, one can develop devices which reflect the specific operating 
parameter of a particular model. A Gas Turbine model for example, could be called 
GT_501K-17 and have all the parameters prespecified for an Allison 501K-17 Gas 


Turbine. 
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4.6 Adding Waveform Types 
Adding a new waveform type requires: 
1. Assignment of a waveform type code. 


2. Writing MATLAB M-Fuile functions for converting to and from the other 


waveform types. 
3, Modification of wconvert .m 


4. Writing MATLAB M-FILE functions for accomplishing the waveform 


operations required by the devices 


oy Modificaiton of wfunction.m files 


a0 








4.6.1 Conversion M-Files 


Here is an example of a conversion M-File for converting a Legendre Series into a 


Polynomial: 

function [poly, jacob] = leg poly (leg,n) 
% 

Syipoly, jacob] = leg poly (leg;,n) 


% Norbert H. Doerry 


% Revision 1.1 21 November 1990 

% 

PmoBourPOlLs sconverts a Legendre Serres to @ Normal Polynomial 
% leg = vector of Legendre Series Coefficients in ascending 
% order 

% n = size of polynomial array to create 
% 

% poly = answer 

% jacob = partial derivative of poly wrt leg 
% 

nl = size(leg); 

mic} = [17 

% 

if n <= 0 

n2 = nl; 
else 

n2 = n; 

end 

% 

% build the jacobian 

% 

jacob = zeros(n2,nl); 

% 

me nl > nz 

nn = n2; 

else 

nn = nl; 

end 


for i=] 2:nn 
jacob(1l:i,1i) = legendre(i-1); 
end 


Oly = jacob * leg; 


dP 3 cP'TD oP 


yee 





Note that a jacobian matrix must be specified for each result with respect to each 
argument. Here is the current version of wconvert .m which is the normal method for 


accessing the conversion routines: 


Bence ton {[w2, Jacob] = wconvert (wl,n,s1,s2) 
% WCONVERT 


VERSION 1.3 of 26 March 1991 
(ere Copyright 1990.1991 by Norbert H. Doerry 


{w2, jacob] = wconvert (wl,n,sl,s2) 


WCONVERT converts a waveform of one type to another type and 
also returns the jacobian of the conversion: 


dP JP dP dP dP dP dP dP dP dP dP dP dP dP oP dP dP dP dP dP dP dP oP oP 


Ww = input waveform 
n = number of points in output waveform 
Si = type of input waveform 
1 data points 
=2 ease series 
= 35 ed bee series 
= 4 ynomial 
= 5 Poe matlab polynomial (not implemented yet) 
= 6 chebyshev series 
$2 = type of output waveform 
nl = size(wl); 
m2) = (); 
meen < 1 
m2 — ni; 
else 
m2 = n; 
end 


% 

meecl < 1 | sli > 6 
waveform type’ 
Ss 
moc Urn) ; 

end 

merseee<. 1 | s2 > 6 
meee waveform type’ 
cS 


return; 
end 
if sl == 5 | s2 == 
is ae waveform type’ 
s 
return; 
end 


Ee 





s2 == 
[w2, jacob] = 
return; 


elseif s2 == 


[w2, jacob] = 


ret urn; 
elseif s2 == 3 

[w2, jacob] = 

TeELuUrNn; 


elseif s2 == 


[w2, Jacob] = 
recurn; 


else 


en 


[w2, jacob] = 
return; 


% 
elceif sl == 


hd 


S == 
[w2, jacob] = 
return; 


elseif s2 == 


[w2, jacob] = 
Setcurn; 
elseif s2 == 3 
[w2, jacob] = 
return; 
elseif s2 == 4 
[w2, jacob] = 
return; 
else 
[w3, 33] = 


w2, JJ 
eb = jj) * 
return; 


end 
% 
elseif sl == 
if s2 = 1 
[w2, jacob] = 
recurn; 
elseif s2 == 2 
[w2, jacob] = 
return; 
elseif s2 == 3 
[w2, jacob] = 
return; 
elseif s2 == 4 
[w2, jacob] = 
return; 
else 
[w2, Jacob] = 
return; 


end 


if s2 == 
[w2, jacob] = 
return; 
elseif s2 == 2 
[w2, jacob] = 
Trecurn; 
elseif s2 == 3 
[w2, jacob] = 
return; 
elseif s2 == 4 
[w2, jacob] = 
return; 
else 
[w2, jacob] = 
return; 
end 


datadata (wl1,n2); 
datafour (wl,n2); 
data leg(wl,n2); 
datapoly (wl1,n2); 


datacheb(wl,n2); 


fourdata (w1,n2); 
fourfour(wi,nZz) 7 
four _leg(wl1,n2) ; 
fourpoly (wl,nZz); 


Lourpoly (wilon2): 
Dae As (woi,71Z)s 
2; 


leg data(wl,n2); 
leg four (wl,nz); 
leg leg(wl1,n2); 

leg_poly(w1,n2); 


leg cheb(wl1,n2); 


polydata(wl,n2); 
polvirour (wl1,1Z) ; 
poly_leg(w1,n2); 
polypoly (wl,n2); 


polycheb (wl1,n2) ; 
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[w2, jacob] = chebdata(wl,n2); 
Tecurn; 

reas ae Z ee 
w3, 33 = chebpoly (wl,0); 
(62, 35] _. = polyfour (w3,n2) ; 
jacob = 35 * 33; 
EreCuUrTh; 

elseif s2 == 
[w2, jacob] = cheb leg(wl,n2); 
return; -_ 

elseif s2 == 4 
[w2, jacob] = chebpoly(wl,n2); 
return: 

else 
[w2, jacob) = leg leg(wl,n2); 
return; = 

end 


end 


Oy GV es 





4.6.2 Waveform Functions 


Waveform functions are defined in a similar manner to the conversion files. Here 1s 


an example of a MATLAB M-File for integrating a polynomial: 
function [p2 , jacob] = poly int (pl,n,c) 


e20, 2 geCOD) = poly anttplpn,;c) 


Norbert H. anes 
Revision 1. December 1990 


input polynomial 
number of points in p2 
integration constant 


pl 
n 
eS 


p2 Imtegeral of pl 
jacob = partial of p2 with respect to pl 


POU INT antegrates a Gaven polynomial and converts the 
result to another polynomial of sizen 


JP dP dP dP dP dP dP oP dP oP oP dP dP OP dP dP OP dP 


oe 


nl = size(pl); 
a (2) = : 


: Peror Checking 


calculate the integration matrix 
ite=— 2eros (nitl,nl); 
or ia = 1 nl 
meeatl, iy = lt) eae 
§ 
: xx is (-1)%*(n-1) 
xX apes (tr BS), 


for i =. J 6 
xx(1,1) = - xx(1,i-1); 


calculate the indefinite integral of the polynomial 
ie 71 * pl; 
convert to a polynomial of the right size 


p27 32] = POLYoOly (pil, nz); 
* ail 


GALI. dP JP DPT 0 dP dO 


% convert to a definite integral by adding the constant of 

% InEeduatton and Subtracting Ehe vaiue of the “indefinite 
: integral at x = -1 
Je 
p2 


Noa eee Ze rOS (2 — ad) | 5 
(1, 1) Oe ee oe a ee Ol, oe ee.) ; 


Ease 





Normally, a user would call w_int to integrate a waveform: 


Lumnceion i[w2 , jacob] = w ant (wl,n,s1,c) 
W_INT 


WMersion is 2 of 19 Apral 1991 
(oe) Copyright. 1991 by Norbert H.. Deerry 


[w2, jacob) = w_int (wl,n,sl,c) 


Weel amtegrares a waveform anG returns the result anto a wavetorm of 
“the same type but of possible different length 


dP de dP dP dP de dP dP de dP dP dP dP de dP de dP de dP dP dP de de dP de 


wl = input waveform 
n = number of points in output waveform 
sl = type of waveform 
= 1 data points 
= 2 fourier series 
= 3 re bees series 
= 4 polynomial 
= 5 for matlab polynomial 
= 6 for chebyshev series 
c¢ = constant of integration 
w2 = waveform which is integral of wl 
jacob = jacobian of w2 with respect to wl 
nl = size(wl); 
mltZ) =. ()7 
iene <1. 
ne = ni; 
else 
n2 = n; 
end 


% 
% check for illegal waveform type 


Pie st <1 | sl > 6 
eartegat waveform type’ 
cS 


return; 
end 
if sl == 
[w2 , jacob] = data int (wl,n2,c); 
elseif sl == 
[w2 , jacob] = four aint (wl,nz,c); 
*lesif $27 Beg poty ut,nd 
w 5 = leg poly (wl,nl) ; 
[w4, 134 = poly int (w3,n2,c); 
[w2, j2] = poly_leg(w4,n2) ; 
HeecOn = 42" 34°": 43; 
elseif sl == 4 
[w2 , jacob] = poly int (wl,n2,c); 
ea. = == 5 ; 7 , 
w 5 = mplypoly (wl,nl); 
[w4, 434 = poly int (w3,n2,c); 
[w2, JZ) = em ee Yate ce) 
FacoD = 32 * 44 * 943; 
elseif sl == 6 
ae 33] = chebpoly(wl1,nl) ; 
w4, 34) = poly int (w3,n2,c); 
[w2, j2] = Boer Wha gs 
HeCob =: 42° ~ 44 -* 43; 
else 
PError’ 
en 


BAG 





Notice how the waveform conversion routines are used to implement waveform 


operations which have not yet been defined for a given waveform type. 


aay 





Chapter 5 Simulation Results 


As a demonstration of the capabilities of WAVESIM, the results of three simulations 
are presented here. While these simulations are relatively simple, they include the important 
features of more difficult simulations, yet are not so complicated as to be unverifiable. The 
first simulation of a simple electical circuit containing only linear devices verifies the ability 
of WAVESIM to construct a viable system and limit truncation error by controlling the 
waveform interval and number of coefficients. The second simulation increases the 
complexity by including a nonlinear device and provides a good test of the Newton-Raphson 
solver. The third and final simulation demonstrates the use of a continuation parameter to 


improve the region of convergence of the simulation. 
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5.1 Linear Electrical Circuit 


To demonstrate WAVESIM’s ability to solve linear circuit problems, the circuit 
shown in figure 5.1-1 was simulated. Initially, both capacitors have zero charge and the 
inductor currents are zero as well. The transients of the capacitor voltages and current are 
shown in figure 5.1-2. The simulation was conducted using Legendre Series with the 20 


second simulation time split up among 23 intervals. Eleven intervals were rejected due to 


excessive truncation error. 


Figure 5.1-1: Linear Electrical Circuit: Schematic 


) R1=.1 ©) Re=l ©) 


© Ce=) 


The results shown in figure 5.1-2 are identical (to working precision) to an analytic 


solution of the circuit. 


Figure 5.1-2: Linear Electrical Circuit: Simulation Results 


rerc.m using LEGENDRE SERIES n =7 


System Variables 
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The input file specifying the system ts given by: 


% 

Pererc. in 

5 

Gevice VDC SOURCE 
TERMINAL 1 
TERMINAL Z 
PARAMETER VS 
END 


% 

device RESISTOR 
TERMINAL 
TERMINAL 
PARAMETER 
END 


WNr 


% 

device RESISTOR 
TERMINAL 
TERMINAL 
PARAMETER 
END 


MNr 


% 

device INDUCTOR 
TERMINAL 
TERMINAL 
PARAMETER 
END 


rFNMoOr 


% 

device INDUCTOR 
TERMINAL 
TERMINAL 
PARAMETER 
END 


Cr Ne 


% 

device CAPACITOR 
TERMINAL 1 
TERMINAL 2 
PARAMETER C 
END 


% 

device CAPACITOR 
TERMINAL 1 
TERMINAL 2 
PARAMETER C 


END 

% 

% 

default 
Gmin 0 
Rmin 0 
EFimport NO 
check both 


Input File for Linear Electrical Circuit 


vo 

1 
0 
0) 


Ral 


ON FF 


R2 


mM Wh 


1b at 


OW 


L2 


oo W 


Ci 


Oon 


C2 


OO W 
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error eqn kel 5e-3 
error eqn pot 5e-3 
error var node  5e-3 
error var flow 5e-3 
Grror malt kelis 10.0 
error mult. pot. “10.0 
error mult node 10.0 
error mult flow 10.0 
max count 10 

max int count 6 

alpha anit. 1.0 

alpha ance anit. .25 
alpha anc min .05 
alpha inc max  .50 
diverge start 3 
diverge max cnt 2 
diverge error mult 10.0 
waveform content max .001 
waveform content nbr 2 
range max .005 

scale potential 1.0 
scale flow IES: 
stype 3 

nor coet 7 

nbr coef min 6 

nbr coef max 14 

nbr data 20 


END 

% 

% 

time 
de min .0) 
dt max 2.0 
Ge Opt «250 
Gt. 2ni1t ~.5 
dt ave 0.0 
start O26 
finish 20. 0 
END 

% 

mot 
potential R1 2 
node 3 
Blow CZ 1 
END 


Figure 5.1-3 shows the time increment and number of coefficients used for each of the 
intervals. Once the transients start to decay, the number of coefficients are decreased to the 


minimum allowed. 
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Figure 5.1-3 Truncation Error Control 


Interval Ending Time (sec) Number of 
Coefficients 


] O25 7 
2 0.50 a 
3 0.75 vi 
4 1.25 7 
S. 225 7 
6 g225 6 
7 4.25 6 
8 4.75 6 
9 525 6 
10 6.25 6 
1] (Ps 6 
12 8.25 6 
13 8.75 6 
14 9.25 6 
i) 10.25 6 
16 M23 6 
17 Epa ae 6 
18 14.25 6 
19 P5225 6 
20 16.25 6 
21 17.25 6 
paps 18.25 6 
23 20.00 6 
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5.2 Nonlinear Electrical Circuit 


To demonstrate WAVESIM’s ability to solve nonlinear circuit problems, the circuit 
shown in figure 5.2-1 was simulated. Initially, the inductor current is zero. As the inductor 
current builds up, its voltage is clamped by the diode to one diode drop above | volt and its 
current ramps up almost linearly. When the inductor voltage falls far enough to turn the 
diode off, the current and voltage both show a normal exponential transient behavior. 


Figure 5.2-2 shows the inductor voltage and current as a function of time. 


Figure 5.2-1 Nonlinear Electrical Circuit: Schematic 


q Ri=1 EQ Di @ 





The results shown 1in figure 5.2-2 were calculated using Legendre Series over seven 
time intervals. Five additional intervals were rejected due to excessive truncation error. 


These results match closely an analytic solution to the circuit. 
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Figure 5.2-2 Nonlinear Electrical Circuit: Simulation Results 


rd.m using LEGENDRE SERIES n = 7 
10 


Pai 


System Variables 
wn 


The input file specifying the system 1s given by: 
Input File for Nonlinear Electrical Circuit 


% 

% rd.in 

% 

@device VDC SOURCE Vsl 
TERMINAL Zt as 
TERMINAL 2 0 
PARAMETER VS 10.0 


END 

% 

device RESISTOR Rl 
TERMINAL - As 
TERMINAL Z 2 
PARAMETER R 1.0 
END 

% 

device INDUCTOR il 
TERMINAL L 2 
TERMINAL Z 0 
PARAMETER L 1.0 
END 

% 

device DIODE1 Di 
TERMINAL 1 Z 
TERMINAL Z 3 
END 

% 
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device VDC SOURCE Vs2 
TERMINAL 1 3 
TERMINAL 2 0 
PARAMETER VS 1.0 
END 


% 

default 
Gmin 0 
Rmin 0 
rimport NO 
check both 
error eqn kcl 5e=35 
error eqn pot Se>3 
error var node 5e-3 
error var flow 5e-3 
error mult kel 10.0 
Srrorsmult pot 10.0 
error mult node 10.0 
error mult flow 10.0 
max count 10 
max int count 6 
alpha init 1.0 
alpha. ane inat 225 
alpha inc min .05 
alpha inc max .50 
diverge start 3 
diverge max cnt 2 
diverge error mult 10.0 
waveform content max .001 
waveform content nbr 2 
range max .005 
scale potential 1.0 
scale flow re 8, 
wtype 3 
Hbr coef -7 
nbr coef min 7 
nbr coef max 10 
nbr data 20 
END 


time 

ae amin 4, 01 
ot max 5.0 
at. Opt... 250 
at) Inte > 
dt ave .0O1 
Start O20 
fanish 10.0 
END 


% 

PICE 
potential Rl 2 
Sow aa) 7: 
END 
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Figure 5.2-3 Truncation Error Control 


Ending Time (sec) Number of 
Coefficients 
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5.3 Nonlinear Mechanical System 


To demonstrate the ability of WAVESIM to use continuation parameters in simulating 
nonlinear mechanical systems, a mechanical power train was modelled. The acceleration 
characteristic of a ship was determined for a propeller rotating at a constant speed. Figure 


5.3-1 shows a schematic diagram of the system. 


Figure 5.3-1 Nonlinear Mechanical System: Schematic 





The propeller model is described in Appendix F-8 while the ship dynamics model is 
described in Appendix F-9. Figure 5.3-2 shows the parametric curves used for C,() and 
C(). This data is for a three bladed propeller with an expanded area ratio of .5 and an H/D 
ratio of .6 [81]. 


Figure 5.3-3 shows the Residual drag coefficient used in the ship dynamics model. 
This data is from the Taylor Standard Series for a hull with beam to draft ratio of 3.0, 
Prismatic Coefficient (C,) of .68, and Volumetric coefficient of 0.002. The Frictional Drag 
Coeffient was calculated using the standard ITTC Line: 


Z TS 
(log io( R.) 7 ay 


CAR.) 
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Figure 5.3-2 Nonlinear Mechanical System: Propeller Characteristics 


3 Blade Propeller Curves: EAR = 5: H/D= 6 


10*CQ and CT 





Figure 5.3-3 Nonlinear Mechanical System: Drag Coefficient 


x10° Residual Drag Coefficient vs Fr 


Cr 





15 = 0.5 0 0.5 1.5 
Froude Number (Fr) 


ral 8.- 





The input file specifying the system is given by: 


Input File for Nonlinear Mechanical System 


% 
= sip mol. an 
% 


& Model of a prime mover attached to prop going to the sea 


% 

device NODE REF GT1 
TERMINAL 1 1 
PARAMETER Vref 5.0 
END 

% 

device PROP1 prop 
TERMINAL SHAFT 1 
TERMINAL WATER 2 
PARAMETER D_ 10.0 
PARAMETER w 0.0 

end 


% 

& parameters are rough ones from Sue B Gail 

% 

device SHIPDYNI1 ship 
TERMINAL WATER 2 
PARAMETER L 100 
PARAMETER A 3300 
PARAMETER M 15000000 
PARAMETER Madd 1.05 
PARAMETER Ca 0.0004 
STATE Us 0 


end 
% 
Node 1 
SCALE POTENTIAL 1 
SCALE FLOW le-4 
NAMES Propeller Shaft 
Gmin 1 
END 
Node 2 
SCALE POTENTIAL 1 
SCALE FLOW le-4 
NAME Hydrodynamic U_ force 
END 
% 
default 
Gmin 0 
Rmin OQ 


rimport NO 

eneck eqn 

error eqn kcl Le-Zz 
CEror, Can peu le-2 
error var node le-2 
error var flow 1le-2 
error mult. kel -10:70 
error multe poet. 10.0 
error mult node 10.0 
error mult flow 10.0 
max count 10 


= he 





Man int count 6 

apna, anit, 075 
alphavine anti 25 
alpha ane min) -—.05 
alpha inc max .50 
diverge start 3 

diverge max cnt 2 
diverge error mult 10.0 
waveform content max .005 
waveform content nbr 2 
range max .005 

scale potential 1.0 
scale flow 10 
stype 3 

nbr coef 7 

nbr coef min 6 

nbr coef max 10 

nbr data 20 

END 


time 

ato many 225 
ae max 5.0 
ac Ope 1.0 
at sine 2:. 0 
at ave 0.0 
Start QO. 0 
Fi naeshs 2/525 
END 


% 
DLOe 
® node 1 converted to RPM 
®%® node 2 converted to Knots (more or less) 
% 
node 1 9.5492966 
node 2 1.8 
flow ship WATER le-5 
flow prop SHAFT le-5 
END 


The results of the simulation using Legendre Series are shown in figure 5.3-4. The 
simulation was broken into 24 intervals shown in figure 5.3-5. An additional 25 intervals 
were rejected due to excessive truncation error. For each iteration, the continuation 
parameter was initially set to 0.5. This value helped assure the initial value of 0 was within 
the convergence region of the nonlinear blocks. While 0.5 was suitable for most iterations, 


several required the continuation parameter be decremented further to achieve convergence. 
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Figure 5.3-4 Nonlinear Mechanical System: Simulation Results 


shipmo2.m using LEGENDRE SERIES n = 7 
50 


45} eat 





35 


30 


System Vanables 
4 


20 knots 


Pe ee eee 
«ee-- 


NY Drag * le-5 N 
5 
O e 
0 2 3 4 5 6 7 8 
time 


As expected, the force on the propeller is greatest during the acceleration of the ship. 
As the ship accelerates, the increased forward velocity on the ship results in smaller torques 
and forces. In reality it is doubtful the motor would be capable of maintaining a constant 


RPM during the acceleration phase. 
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Figure 5.3-5 Truncation Error Control 


Interval Ending Time (sec) Number of 
Coefficients 














l ] 

2 i Re 8 
3 1-75 10 
6 Pap bs 10 
> 25 10 
6 PRES. 10 
i 3.0 10 
8 2) 10 
9 oe) 10 
10 a5 10 
Jet 4.0 10 
| Ps 4.25 10 
13 4.5 10 
14 4.75 10 
1S 5.0 10 





16 5,25 10 
17 SB. 10 
18 TEU) 10 
19 6.0 10 
20 6.25 10 
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—_" 
a 
WN 
— 
_) 





- 162 - 





Chapter 6 Conclusions 


In its present form, WAVESIM is ideally suited for testing numerical algorithms. 
While it is capable of simulating large systems, the interpretive nature of MATLAB is not 
numerically efficient enough for serious simulations. Careful development of a simulation 
environment based on the techniques explored in WAVESIM should prove effective in 


solving tightly coupled multirate systems of lumped parameter models. 


The simulation environment described in this thesis should be considered a framework 
for future developments. Many improvements are possible and desirable. In particular, the 


following areas need further attention: 
Truncation Error Control 


The present method for controlling truncation error is heuristic and should be 
examined for improvement. Truncation Error propagation should be examined and 


given a theoretical basis. 
Discontinuity Time Prediction 


The accuracy of the methods used in WAVESIM depend partly on the ability to 
predict discontinuities and force them to occur on time interval boundaries. The 
methods used in current models are crude and should be replaced with more robust and 


accurate methods. 
Stability Analysis 
WAVESIM presently does not perform any stability analysis. Since WAVESIM 


abandons the standard state space representation of the system, determining the 
eigenstructure of the system is not easy. A Stability measure based on the 
characteristics of individual devices would fit well with the structure of WAVESIM 
and would be quite useful in the design of distributed controls. 


Smoothing Operation 


The smoothing operator for removing the effects of high frequency 
discontinuities needs to be examined to improve its efficiency. How long to make the 


smoothing interval is a question which has not been satisfactorily answered. 


Partitioning and Relaxation 
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The approach WAVESIM uses for developing the set of system equations is 
ideally suited for use with relaxation methods if the system is weakly coupled. An 
extension to the structural Jacobian to include matrix norms would greatly simply the 
task of partitioning the system into a set of weakly coupled blocks which internally are 
strongly coupled. Each individual block would be solved using Newton-Raphson with 
the system solved using a relaxation technique. Unfortunately, the process of 
constructing a system matrix of norms from device matrix norms is presently not 


possible because arithmetic operators for matrix norms have not been identified. 


Overall, WAVESIM has been very successful in developing the algorithms for building 
Systems in terms of device functions, treating waveforms as an abstract data type, and 
employing the structural Jacobian matrix to reduce the system into a sequence of smaller 
blocks. Much work remains, but the foundation of a waveform based simulator capable of 


handling tightly coupled multirate simulation problems 1s contained within WAVESIM. 


PGA 
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Appendix A: Glossary 


Block 


Continuation 
Parameter 


Device 


device.def 


Device Jacobian 


Device Structural 
Jacobian Matrix 


Export Variable 


A subset of a system’s equations and variables which must be solved 
simultaneously. Blocks are organized into a sequence where 
variables determined in a previous block may be used in following 
blocks. 


A technique to enlarge the convergence region of a nonlinear system 
by using the solution to a linear system as the initial guess for the 
solution of another system which is a combination of the linear and 
nonlinear systems. The process is repeated with each iteration 
increasing the nonlinear portion until the solution to the nonlinear 
system is determined. The continuation parameter determines the 
relative proportion of the nonlinear system: O for the linear system 
and 1 for the desired nonlinear system. 


A device is a mathematical model of a physical piece of equipment 
comprising a system. Devices interact with one another through 
interface variables which are associated with other device interface 
variables through terminals connected at nodes. The equations 
describing a given device type are specified in the device definition. 
A given instance of a device also has associated parameters and 
nodal connections. 


A file for descnbing a device definition. Each device type has an 
entry describing the device type name, terminals, states, parameters, 
structural Jacobian, and MATLAB M-Fuie containing the 
constitutive equations. 


A matrix whose elements are the partial derivatives of the export 
variables of a device with respect to the device import variables. 


A matrix describing the dependence of a device’s export variable 
with respect to the import variables. The dependence is specified by 
a matrix whose elements are a code indicating if the dependence is 
zero, identity, diagonal, linear, or nonlinear. 


An interface variable (either a potential or flow variable) of a device 
which is explicitly defined by the device constitutive equations. A 
device takes import variables as input and produces export variables. 
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Fiow Variable 


ca 


Import Variable 


Information 
Terminal 


Interface Variable 


KCL 

KCL Equation 

KCL Group 
Number 


MATLAB M-File 


Newton-Raphson 
Method 


An interface variable (either an export or import variable) associated 
with a normal terminal of a device which corresponds to a quantity 
satisfying Kirchhoff’s Current Law at nodes. Examples of flow 
variables are currents, forces, and torques. 


A modification to the KCL equations at a node corresponding to the 
insertion of a conductance to the Q potential. Used to prevent 
singular systems. 


An interface variable (either a potential or flow variable) of a device 
which 1s implicitly defined by the device constitutive equations. A 
device takes import variables as input and produces export variables. 


A terminal of a device having only a potential associated with it. 
Used to convey energyless information between devices. 


Variables through which devices communicate energy and 
information transfer to other devices. Interface variables are 
associated with terminals, can be classified as either flow or 
potential variables and can be classified as either import or export 
variables. 


Kirchhoff’s Current Law which states the sum of the flow variables 
attached to a node is zero. 


Equates the sum of the flow variables attached to a node to zero. 


If a subset of a device’s flow variables add to zero by definition, 
then the elements of such a subset have a device-unique nonzero 
group number. Flow variables which do not belong to such a subset 
have a 0 KCL Group Number. 

KCL Group Numbers are used to determine possible singular 
systems due to linear dependence of system KCL equations. 


Text files of MATLAB commands for creating new MATLAB 
functions or executing Scripts. 


An iterative technique for solving systems of nonlinear equations 
which uses the Jacobian Matrix to generate corrections to the system 
variables. 
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Node 


Normal Terminal 


Parameter 


Potential 
Difference 
Equation 


Potential Variable 


fe, 


smoothing 
Operator 


State 


A connection point for connecting terminals of one or more devices. 
If at least one normal terminal is attached, the node is a normal node 
and a system KCL equation is written to equate the sum of all the 
attached flow variables to zero. If only information terminals are 
attached to a node, the node is an information node. All nodes have 
an associated node potential. 


A terminal having both a potential and flow variable. Used to 


simulate energy transfer between devices. 


A variable which does not change throughout the simulation. 
Usually refers to machine ratings, resistances, time constants, etc. 


Each export potential variable in a system has an associated 
potential difference equation equating to zero the difference between 
the potential of the node to which the variable is attached and the 
value of the export potential. 


An interface variable (either an import or export variable) associated 
with either a normal or information node. Al of the potential 
variables attached to the same node are equal to the potential of the 
node. 


A modification to a potential difference equation corresponding to 
the insertion of a series resistance. Used to prevent singular 
systems. 


A waveform operator for removing the high order waveform content 
of its argument by retuming a waveform which is the convolution of 
the argument with a square pulse. The returned waveform is 
effectively the local average of the argument waveform. The 
smoothing operator removes unnecessary detail from a waveform 
and improves the representation of the desired properties of the 
waveform with fewer coefficients. 


A State 1s a device variable whose value is retained between 
waveform intervals. Constants of integration and device operating 
modes are the most common uses of states. 
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Structural Jacobian A code indicating the nature of an element of a Jacobian matrix: 


Code 


Subsystem 


System 


System Jacobian 


System Structural 
Jacobian Matrix 


System Variable 


Terminal 


Waveform 


Y/aveform Content 


Zero Matrix 

Identity Matnx 
Diagonal Matrix 

Linear Matrix 

Type A nonlinear Matrix 
Nonlinear Matrix 
Unknown 


CA en eee 


A subset of the devices of a system which are grouped together and 
solved independently of other devices and subsystems. Subsystems 
have not been implemented in WAVESIM. 


A group of devices and subsystems and the nodes interconnected 
them. 


A matrix containing the partial derivatives of the system equations 
with respect to the system variables. 


A matrix describing the dependence of a system’s equations with 
respect to the system variables. The dependence is specified by a 
matrix whose elements are a code indicating if the dependence is 
zero, identity, diagonal, linear, or nonlinear. 


The set of system variables 1s composed of all the node potentials 
and all the device import flow variables. 


A modelling analogy to a physical attachment point on a device. 
Normal terminals have an associated flow and potential variable and 
are used to model the transfer of energy into and out of a device. 
Information terminals have only a potential variable and are used to 
convey information between devices. 


A representation of a variable over a given time interval consisting 
of a vector of coefficients and a waveform type indicator for 
specifying how the coefficients should be interpreted. Common 
waveform types are Legendre Series, Chebyshev Series, 
Polynomials and Data Points. 


The magnitude of a coefficient of a waveform divided by the square 
root of the sum of all the waveform coefficients. The Waveform 
Contents of the higher order coefficients are used to determine if the 
truncation error is negligible. 
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Waveform Type An indicator specifying how the coefficients of the waveform vector 
Indicator should be interpreted. Common waveform types are Legendre 
Series, Chebyshev Series, Polynomials and Data Points. 


WAVESIM A numerical algorithm development program incorporating the 
systematic treatment of waveforms as a data type, the terminal 
description of devices, and the use of structural Jacobians in system 
reduction. 


gg 








Appendix B: Continuation Parameter Pitfalls 


If used properly, continuation parameters can help enlargen the region of convergence 


of an iterative scheme. This section will show how continuation parameters can fail due to 
bifurcations of solutions. 

Take for example, the following system of two equations and two unknowns 
Ftxy) = 0: 


F(x, y)=y —@° —x)=0 
F(x, y)= y —(Mx+B)=0 


Initially, set M = 0 and B = 1.875. From the following figure, it 1s obvious the solution 


is the intersection of the two curves and falls at the point (1.5,1.875). 


Figure B-1: Solution to y = x° - x and y = 1.875 


ae Cn SS ee 


Nonlinear Solution 





To solve this system with a continuation parameter, we create a new function H(x,y,q) 


which is formed by combining F(x,y) with a linear system G(x,y): 


G(x, y)=y —(mx +b) = 0 


G,(x,y) = y —(Mx +B) =0 
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Hy, O) = OF Coy) (1 O)Gx; y)=0 


The modeller now has the choice of selecting m and 6. A natural choice would be a 
linearization about a given point. If we linearize about x = 0, the values are m = -1 and b = 0. 


The following figure shows the results of this selection: 


Figure B-2: Continuation Method for m=-1 b=0 M=0 B=1.875 


y = a(x°-x)+(1-a)(-x) 





ROOT LOCUS for y = 1.875 y = a(x43-x) + (a-1}(-x) 
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a = 


x Root Locus Points 







y =(x°-x)+(1 — a) (-x) 
ES 


ec 


Note the solution for @=0 is (-1.875,1.875) which is not very close to the desired 
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solution for @=1. Furthermore, as © increases slightly, it actually becomes slightly more 
negative until the nonlinear curve no longer intersects the linear equation in the left hand 
plane. At this point, the solution has a discontinuity and jumps into the right hand plane with 
a value for x much larger than the solution. The root locus for x as a goes from 0 to 1 clearly 
shows this. Hence for this selection of m and BD, the use of the continuation parameter makes 


the job of solving the system tougher instead of easier. 


If we choose different values for m and Bb, the situtation may change. Say for example, 
we set m=I1 and b=90. This selection appears to work well as can be seen with the 


following figure: 
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Figure B-3: Continuation Method for m=1 b=0 M=0 B=1.875 





y = a(x?-x)+(1-a)x 





ROOT LOCUS for y = 1.875 y = a(x43-x) + (a-1}){x) 
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x Root Locus Points 


y= o(x° —x)+(1-a)x 


y = 1.875 


1.7891 -0.8945 - 5.03991 -0.8945 + 5.03991 
1.6440 -0.8220 - 2.24211 -0.8220 + 2.24211 | 


-0.7768 - 1.34551 -0.7768 + 1. pee 





170.7500 - 0.82921 -0.7500 + 0.8292 


The solution for @=0 is close to the solution and as & increases, it rapidly converges 
on the desired solution (1.5,1.875). We should not rejoice however, because even this 
selection can fail for other choices for M and B. For example, if M = 2.5 and B = -3, the 
following figure demonstrates a discontinuity in the solution path: 


Figure B-4: Continuation Method for m=1 b=0 M=2.5 B=-3 











= 152. 








ROOT LOCUS for y = 2.5x-3 y = a(x43-x) + (a-1)(x) 


x Root Locus Points 


y =a(x°—x)+(1-a)x 
y=2.5y-3 


a 


1.3445 - 0.65071 1.3445 + 0.65071 


0.0735 -5.4659 2.7706 2.6953 


1.1024 - 0.38151 1.1024 + 0.38151 





1.0000 
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The bottom line is that it may not be possible to develop a transformation function 
whose solution vector is always continuous. Any information known as to the region where 
the probable operating point lies should be used in directing the solution to that region. For 
this example, if x is known to be constrained to the interval [-5 5] and M 1s known to be less 
than 17.75 (.75¢5° - 1 is the slope of the line tangent to y = x° - x and passing through (5,120) 
) then y is also constrained to the interval [-120 120]. If we use as our linearizing function 
the line connecting (-5,-120) and (5,120) it is clear the root locus will also remain within the 


constraints for any value of M or B meeting the constraints at a= 1: 
y =x —x)+(1-0)24x 


Figure B-5: Continuatin method for m=24 b=0 A1=2.5 B=-3 


y = a(x?-x)+(1-a)24x 
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ROOT LOCUS for y = 2.5x-3 y = a(x43-x) + (a-1)(24x) 





x Root Locus Points 
y =ax>—x)+(1—0)24x 
y=2.5y-3 


0.9500 0.9945 + 0.77371 


1.0000 1.1024 - 0.38151 1.1024 + 0.38151 


0.5000 0.1657 - 4.25231 0.1657 + 4.25231 





0.6000 
0.7500 
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For this example it is actually quit easy to determine if a bifurcation will occur. At a 


bifurcation, the x root locus points satisfy the following relationship: 
(xc) (x-d)=0 
x? +(-d —2c)x? + (c? + 2cd)x + (-c’d) 


where c is the multiple root whose paths will deviate from the real x-axis and d is the 


root staying on the x-axis. Now the actual equation defining the roots 1s given by: 
ox? +((1-a)m —a-M)x+(1-a@)b -B =0 
Equating terms we get: 


—d—-2c =0 


1 
c* + 2cd =~ ((1— om —a—M) 


24-4 _ yy 
-c'd =—((1—a)b -B) 


solving this system for c, d, and Q, we get: 


d =-—2c 


mil =O) 2 
a 20 


(l-o)b-BY_ 
(1-—a)m -a-M+3o{ C=G2=8 | =U 


If one of the solutions for @ is a real number in the interval [0,1], then there will be a 
bifurcation and possibly a discontinuity in the path. If two of the roots approach from +c 
and -co along the real axis and a solution for & exists in the interval [0,1], then there will 
definitely be a discontinuity in the solution path. If two roots appraoch from off the real axis, 
combine at the bifurcation point, then travel in the +x and -x directions, there will be three 
real solutions for x and the solution path will converge onto one of them. If there is no real 
solution for & in the interval [0,1], then there will be no bifurcation, no discontinuity in the 


solution path and the solution will be unique. 
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Appendix C: Load Flow Example 


The method for building systems can be applied to static simulations for determining 
equilibrium points of systems. The traditional load flow is representative of this type of 
problem. Figure C-1 shows a three bus load flow example consisting of four device types: 


PV Generator, VD Generator (Slack Bus), PQ Load, and a transmission line. 


Figure C-1: 3 Bus Load Flow Example 
Niles 
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C-1:;: Device Definitions 


C-1.1: PV Generator 


Interface Variables 


Terminal Potential Variable Flow Variable (KCL Group) Type 
VQ V (export) Q (import) (0) Normal 
DP D (import) P (export) (0) Normal 


The import X;,,,, and export X,,, vectors are defined by: 


~[e 





: Parameters 
| 
| ee Scheduled Generator Power 
Vo Scheduled Generator Voltage 
Equations 
V=V. 
j ey ge 


Device Structural Jacobian 


The device structural jacobian is given by: 


0 0 
Inox 


Device Jacobian 


The device jacobian is given by: 
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C-1.2: VD Generator (Slack Bus) 


Interface Variables 


Terminal Potential Variable Flow Variable (KCL Group) Type 
VQ V (export) Q (import) (0) Normal 
DP D (export) P (import) (0) Normal 


The import X;,,, and export X,,, vectors are defined by: 


Le 


V 
Xan =| 5 | 


Parameters 


Ve Scheduled Generator Voltage 
Dz Scheduled Generator Angle 


Equations 


Device Structural Jacobian 


The device structural jacobian is given by: 


0 0 
Jos=| g 


Device Jacobian 


The device jacobian is given by: 


0 0 
‘orto 0 
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C-1.3: PQ Load 


Interface Variables 


Terminal Potential Variable Flow Variable (KCL Group) Type 
VQ V (import) Q (export) (0) Normal 
DP D (import) P (export) (0) Normal 


The import X;,,, and export X,,, vectors are defined by: 


wl 


ne 
ZAlde 


Parameters 


lige Scheduled Load Real Power 
QO, Scheduled Load Reactive Power 


Equations 


Device Structural Jacobian 


The device structural jacobian is given by: 


0 0 
Jos=| 9 


Device Jacobian 


The device jacobian is given by: 


0 0 
ola 0 
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C-1.4: Transmission Line 


Interface Variables 


Terminal 


VQl 
DPI 


vQ2 
DP2 


Potential Variable 


V, GQmport) 
D, (import) 


V, (amport) 
D, (import) 


Flow Variable 


Q, (export) 
P, (export) 


Q, (export) 
P, (export) 


(KCL Group) Type 


(0) Normal 
(0) Normal 


(0) Normal 
(0) Normal 


The import X;,,, and export X,,, vectors are defined by: 


V, 
D, 
V, 


WN 


D 


WN 


exp oF 


Parameters 


R Transmission Line resistance 
ne Transmission Line reactance 


Equations 
Obtain Y: 


~ R 
R24? 


Ems 95 
R? +X? 
Y = VA?+B? 


Dy = atan2(B,A) 
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Calculate Side one current 
= Vecos( Dt D1.) cos(D, + D,) 


I, =V,Y sin, + Dy) -V,Y sin(D, + Dy) 


Calculate Side two current 


Calculate real and imaginary parts of the voltages 
Viz = V,cos(D,) 
V,,=V,sin(,) 
Vip = V,cos(D,) 
V,, = V,sin(D,) 
Calculate the export variables (Powers) 
ey lane gly 
0, aa Vinly 1 Virtir 
P= Von lip ck Var 
O, =—Voplo, + Vazlop 


Device Structural Jacobian 


a ee en ee 
Seo a 
Sa 
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Device Jacobian 


Calculate the Partial derivatives of the voltages with respect to the import variables: 


Fe = eos) -V,sinD,) 0 Oj 





wv, 


ax, Sunt? v V,cos(D,) 0 0] 





WVoq 
ax, 


imp 


OV, 
OX; 


imp 





=(0 0 cos,) —V,sin@,)] 


=10 5°00 2sin(25)4 V5 cos.) 





Calculate the partials of the currents with respect to the import variables: 

















ol 
a =[Ycos(D,+D,y) —YV,sin(D,+Dy) -Ycos(D,+Dy,) YV,sin(D,+Dy)] 
imp 
aly 
ya [Ysin(D,+Dy) YV,cos(D,+Dy) -—Ysin(D,+D;) -—YV,cos(D,+Dy)] 
imp 
np Ap 
OX ON, 
ox Ox, 


imp imp 
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Calculate the jacobian matrix 


Vg 
ax 


imp 
Vy 
sp 
OX imp 


Ce |e) 
avy Foe 0 0 0 ior 
V = 








imp 

Of Sh. ie Ve = Valores 
Te all Roxeee 

avs; 

Oe 

ane 

OXinp 

a 

OX imp 


oy 
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C-2: Network Description 


Figure C-2-1 details the device interconnections of the 3 Bus system shown in Figure 
G1. 


Figure C-2-1: 3 Bus Loadflow Block Diagram 
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C-2.1; Variable Labeling Convention 


For this example, the following convention will be used for labeling variables and 


functions: 
Device Terminal Variables: x,, 55 ca 


aa Device Name 
bb Variable Name 


Cc n = normal terminal 
i = information terminal 


d t= import variable 
e = export variable 


f p = potential variable 
f= flow variable 


Device import variable vector x,, ; 
aa Device name 
Device export variable defining function x44 55 cay = Soa bb caf\Xaa i) 


aa Device Name 
bb Variable Name 


Cc n = normal terminal 
i = information terminal 


= export variable 


Sf p = potential variable 
f= flow variable 


Device Jacobian J,, 


Device Jacobian Element J. 55 ¢¢ 


aa Device Name 

bb ‘Export Variable Name 

gg Import Variable Name 

System Variables: Node Potentials V, 
n Node Serial Number 


System Variables: Flow Variables J,, ;, 


aa Device Name 
bb Variable Name 
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System Equation: KCL g,() 
n Serial number of node KCL 1s applied to 
System Equation: Potentials g, .. »() 


n Serial number of node 
aa Device Name 
bb __ Export Potential Variable Name 


System Jacobian Element: KCL vs Node Potential J,,, , ,, 

System Jacobian Element: KCL vs Import Flow J,,, , oa 5 

System Jacobian Element: Potential Eqn vs Node Potential J,,, .. aa m 
System Jacobian Element: Potential Eqn vs Import flow J,,, .< aa aa bb 
n Serial number of KCL node 

m Serial number of Node Potential 

aa Flow Variable Device Name 

bb Flow Variable Device Variable Name 


GC Potential Equation Potential Device Name 
dd _—~ Potential Equation Potential Variable Name 
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C-2.2: Network Specification 


Now that the variable labeling convention has been addressed, it 1s time to define the 


devices and the network interconnecting them. 


PD Generator Gl 


Terminal Potential Variable Flow Variable Node 
DP fone co 
Parameters 

Ve 1.0 PU 

De 0.0 RAD 

Import Vector: 


= XG1_O nif | I61 0 
Gli = 
AGi_P nif Ie) p 


PV Generator G2 

Terminal Potential Variable Flow Variable Node 
MAS. XG2_V_nep XG2_Q nif 2 

DP XG2_D_nip x G2_P_nef 5 
Parameters 

| ae 0.5 PU 

Ve 1.05 PU 

Import Vector: 


en XG2_Q nif | _ I62 0 
G2_i io 
NGoD nie Vs 
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PQ Load L3 


Terminal Potential Variable 
VQ X13 V_nip 

DP X13 D_nip 
Parameters 

P, 0.6 PU 

0, 0.3 PU 

Import Vector: 


Transmission Line T12 


Terminal Potential Variable 
VQl X712_V1_nip 

DP1 XT12_D1_nip 

VQ2 X112_V2_nip 

DP2 X112_D2_nip 
Parameters 

R ‘em a6) 

4 0.60 PU 

Import Vector: 


ES oo ae 


Flow Variable Node 
X13 _O_nef 3 
X13 P_nef 6 
ey V; 
X13 _D_nip Ve 
Flow Variable Node 
XT12 Q1_nef l 
X1T12 P1_nef 4 
Xr12_02_nef Z 
XT12_P2_nef 2 
XT12_V1_nip Vy 
2a I snip.|| 2 V, 
X7}2_V2_nip Vy 
XT12_D2_nip Vs 
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Transmission Line T13 


Terminal Potential Variable Flow Variable Node 
VQ! X713 V1_nip X713 Q1_nef 1 
DP1 X713_D1_nip X713 PI_nef 4 
VQ2 X713_V2_nip X713_Q2_nef 3 
P2 X713_D2_nip XT13_P2_nef 6 
Parameters 
R 05-20 
xX 0.20 PU 
Import Vector: 
XT13_V1_nip Vi 
XT}3_D1_nip V 
MEI = a 
XT}3_V2_nip V; 
X13 D2 nip Ve 
Transmission Line T23 
Terminal Potential Variable Flow Variable Node 
VQ! X723_V1_nip X723_Q1_nef 2 
DPI X723 D1_nip X7123 P1_nef 5 
VQ2 X723_V2_nip X123 Q2_nef 3 
DP2 X123_D2_nip X7123_P2_nef 6 
Parameters 
R 0.10 PU 
xX 0.40 PU 
Import Vector: 
X723_V1_nip V, 
Pe a X73 pinip| _| Vs 
i X123_V2_nip V; 
X723_D2_nip Ve 
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C-2.3: System Variables and Equations 


There are nine system variable and equations associated with this example. There 
are the six node potentials plus three import flow variables ordered in the following 


manner: 
7 
ee 1A Vue ee ee Ion 0 Io) p loo g] 


The nine system equations are composed of six Kirchhoff Current Law equations and 


three potential equations: 
iC ae) = 16) 9 +X712 91 nef +X713_O1 nef 
£4(X,,.) = 162 9 +712 92 nef +X123_91 nef 
£3(%5) = XL3_9 net + 113_02_nef +%123_02 nef 
e se) =16) p + Xr12 pt_nef + %713_P1_nef 
g 5(X,,,) =XG2_ p nef +X112_p2_nef + *123_P1_nef 
c sen) =X13 p_nef t%113_p2 nef tX123_p2_nef 
§ 1cr_v%gys) =V, Ser GLY ee 
25 G2 v%qys) =V,- XG2_V_nep 


SAnGIED oe) =V,- XG1_D_nep 


- 201 - 





C-2.4: System Structural Jacobian Matrix 
The equations for generating the system Jacobian matrix are given by: 
£1(X,55) 


J ys =4 112 01-v1 +4713 91. v1 =N +N 


« 


Deep =/T12 Q1_V2 = 


7 


J sales =J 713 91-y2 =N 


a 


J ys1.4=4 712 91-p1 t+ 4713.91 pp =N +N 


Jys15 =! r2 501 2S N 


- 


aie = 713 0! Dp? = N 


J ys 1_G1_O = | 
Px.) 
J js2 1 = J 712 502.) — N 


J ys 22 = Pie02 V2 + Jiro; oly > N+N 


J 52.8 = J 793 SONS VO = N 


J ys 24 = 717 62:)1 = N 


J ys.2.5 =4 712 92.2 +4123_91p1 =N +N 


ae =J1>3 91 p? =N 


J fys22G200: i 
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Pa(%,,,) 


J sys_3_1— J 773 SO7Vl = N 
J sys_3.2 —~%123_92_V1 = N 


J s3.3 =I 13 9 v tI 113 92. v2 +473 92.2 =OtN+N 


- 


J oe) tae S743 oa 0 em N 


J sys_3_5 ~%T23_92_D1 — N 


Uineste = 5000 +Jr13 92 p2 +J3 02 p? =O0+N+N 


g& a Xqy5) 


J sys4.1 =! 112 p1_vi +4743 piv) =N +N 


J sys_4_2 — J Tia Pia N 


J Cy te ie J TIS PIV N 


- 


J sys_4_4— s TI2_P1_p1 + J TI3 Pipl = DNerUN. 


J sys_4_5 —~%T12_P1_D2 — N 
J ys4.6 =! 713 p1_p2 = 


sy 


ual Pp a if 
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aX) 
J re ae J TI2P2 VE N 


J svs 52 a T12_P2_v2 t J 73 Pieve = N+N 


J sys_5_3 — J TORII oe N 


J sys_5_4 —4T1I2_P2_D1] — N 


J sss = 462 p_pt4r12_p2_p2 +4123 p1_p1 = ON ate 


- sys_5_6 — %123_P1_D2 — N 


J ys G20 Glo 0 


© 


&AXsys) 
J noe ame 0 ae N 
J sys_6.2 — 4 T23_P2_VI — N 


J ys6.3 = 413 pv tI 113 p2_v2 tI103_p2 v2 =OtN+N 


J sys 6.4 — J Ti35P2 Di = N 


J sys_6_5 — J 723_P2_D1 — N 


J ys66 = 413_P_p +473 p2_p2t4703_ p2 pr =DtN+N 
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£1 cv X45) 
J sys_G1_V_1— I 


J ys_G1V_G1_Q = G1_V_9 — 0 


J os GIV_GI_P = o1_v_p =9 


« 


82 Gz VX qy5) 
J sys_G2_V_2~ is 


J ys G2V_5 = —J G2_V_D = 0 


J ys G2_V_G@2 Oe =I, Gv = 0 


£4 G1 _p(Xsys) 
J os GI _D_4 =[ 


J ys GI_D_G1_Q = Jo: poe = 4) 


J sys_G1_D_G1_P = /G1_p_p =9 


« 
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oo 


By applying the rules of Structural Jacobian element arithmetic on the system 


equations, we can generate the following system Structural Jacobian: 


ah eee ee, ee 
Gt ee eee 
cooo2z22222 
OS ee aes ee ee 
ooo z2z2zZzz aa 
S21) ea ee 
Sia2o. = /Sfeveonar 
S1O6e%S oe oe 1 
Sere © 1 o4eleo =o 


Close inspection of this matrix reveals seven blocks: Six 1X1 element blocks and one 


3x3 element block: 
Block 1 
System Row: i 
System Column: ] 
System Variable: V, 


Equations: 
Biteny te) = V, —~%XE1_V_nep 
Structural Jacobian: 


J, =(] 
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Block 2 
System Row: § 
System Columns: 2 
System Variable: y, 
Equations: 
8 G2_V\%qys) = Vo —XG2_v_nep 
Structural Jacobian: 
Ja | 
Block 3 
System Row: 9 
System Column: 4 
System Variable: Ve 
Equations: 
84 G1 pDXys) = Vo Gi D_nep 
Structural Jacobian: 


Jp3 = [1] 
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Block 4 
System Rows: 3 5 6 
System Columns: 3 5 6 
System Variables: V; Vs, V, 
Equations: 
£3(X,,.) = X13 0 nef + *713_02_ nef + X723_02 nef 


85(%,,5) = %G2_pnep + %112_P2_nef + X123_P)_nef 


K Oe) =Xr3 p nef t%113_p2_nef + X123_P2_nef 


Structural Jacobian: 


Block 5 
System Row: J 
System Column: fi 
System Variable: IG10 


Equations: 
g jee) =16) 9 + X71 91 nef +*713_01 nef 
Structural Jacobian: 


Jes =[] 
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Block 6 


System Row: 2 
System Column: 8 
System Variable: Ler p 
Equations: 


84(%,,5) = 162.9 + X712_92 nef +%123_91 nef 


Structural Jacobian: 


Jp6 = [1] 
Block 7 
System Row: 4 
System Column: ) 
System Variable: I G2 9 
Equations: 


B4(% ys) = 4161p + Xr12_p1 nef +X713_P1_nef 
Structural Jacobian: 


Ja7 =) 
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C-2.5: Solving the System 


Applying the equations for the first three blocks yields: 


0 
Vi = foveal] 4) = 1.0 PU 


V,= for. aa | 4) =1.05 PU 


0 
V, = fa. 0e( |) = 0.0 PU 


Now the following system of three equations for the fourth block must be solved: 


83(Xp4>Xpre) 
Meroe = £5(Xpa>Xpre) - 0 


§ 6(Xp4 ’ ee) 


Where: 
Ve 
Xpq = Vs 
Ve 
Vi 
X ore = V2 
V, 


Starting with the intial guess of [1 0 0]' for x,, we obtain the following error vector 
and jacobian matrix: 


0.1824 
Xror =| — 0.4485 
0.5706 


6.9412 06176 —-1.7941 
Jepg=| —0.6176 4.1176 —2.4706 
1.7253 -2.4706 7.1765 
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Inverting the Jacobian and multiplying by the error vector results in the following 


correction vector for Xp, : 


0.0441 
x, =| —0.0769 
0.0424 


i 0 0 
Xpq = Xpq Xa 


By repeating the Newton-Raphson iterations several more times, the following table 


can be constructed: 


i 1.0000! 0.00001 coool 0.18241 ~~ -0.4485 | 0.5706 


ay 0.9502 0.0762 -0.0463 = 0.448e-4 | 2.208e-4 
31 oss0a} oarei| -o0ssa] 137028] -oniea|  1200re 


From these results, the final three blocks can easily be solved: 












IG) 9 = fri 91 nef%r12_1) —Sry3 01 nefX113 5) 


T6109 = -0.1451 PU 


Loy p = Frir_pr_nef%112.i1) — Sti3_p nef113_i) 


Ig, p = -9.1233 PU 


IG 58 laos Shri _02-neAXr12_i) — frr3 01 _nef\X123_) 


I¢2 g = -9.2483 PU 
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C-3: Summary of Results 


Bus Voltage Magnitude 
Bus Voltage Angle 


Gl Real/Reactive Power 
T12 Real/Reactive Power 
T13 Real/Reactive Power 


Bus Voltage Magnitude 
Bus Voltage Angle 


G2 Real/Reactive Power 
T12 Real/Reactive Power 
T23 Real/Reactive Power 


Bus Voltage Magnitude 
Bus Voltage Angle 


L3 Real/Reactive Power 
T13 Real/Reactive Power 
T23 Real/Reactive Power 


Bus I 


1.0 PU 
0.0 rad 


-0.1233 PU 
-0.1437 PU 
0207 1rd 


Bus 2 


1.0500 PU 
0.0761 rad 


-0.5 PU 
0.1471 PU 
03529 PU 


Bus 3 


0.9502 PU 
-0.0464 rad 


0.6000 PU 
-0.2617 PU 
=()'3333 PU 
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-0.1451 PU 
-0.0423 PU 
0.1874 PU 


-0.2483 PU 
0.0558 PU 
1925, PU 


0.3000 PU 
-Q.1661 PU 
“0.1339 PU 





Appendix D: Modified Load Flow Example 


Appendix C demonstrated how a system can be built and solved for a conventional 
load flow problem. This example demonstrates how control signals such as real and reactive 
power sharing signals can be incorporated in the load flow solution. In particular, this 
example connects two parallel generators to a load via a transmission line. A conventional 
load flow fails for this example because the generator bus voltage magnitude is 
overdetermined and there is no relationship for sharing reactive power. In this example, 
information variables are used to force each generator to be proportionally loaded and have 


the same power angle. 


Figure D-1: Parallel Generator Load Flow Example 


o 
des 
°C 


G2 


D-1: Device Definitions 


In addition to the transmission line and PQ load defined in Appendix C, two more 
devices must be defined: A slack bus generator incorporating the load sharing information, 


and a PQ generator employing the load sharing. 


= 1S 





D-1.1: VDS Generator (Slack Bus) 


Interface Variables 


Terminal Potential Variable Flow Variable (KCL Group) Type 
VQ V (export) Q (import) (0) Normal 

DP D (export) P (import) (0) Normal 

Pp p (export) Information 

q g (export) Information 


The import x;,,,. and export x,,, vectors are defined by: 


ae u 
ep P 1D 
Xap = ' 
q 
Parameters 
Ve Scheduled Generator Voltage 
De Scheduled Generator Angle 
P, Scheduled Generator Power Base 
Equations 
V=V, 
D=D, 
oe 
_Q 
a. 


Device Structural Jacobian 


The device structural jacobian is given by: 


Jos = 


Pe aye = 
2 o-oo eS 


aA 





Device Jacobian 


The device jacobian is given by: 


0 Q 
Q QO 
1 
oe 
lp go 
P ~ pil 
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D-1.2: PQS Generator 


Interface Variables 


Terminal Potential Variable Flow Variable (KCL Group) Type 
VQ V (import) Q (export) (0) Normal 

DP D (import) P (export) (0) Normal 

p p (import) Information 

q q (import) Information 


The import x;,,,, and export x,,, vectors are defined by: 


V ee QO 
D 2a 2 
Nini = 
se 
q 
Parameters 
jeg Scheduled Generator Power Base 
Equations 
P=-P;p 
QO =-P,pq 


Device Structural Jacobian 


The device structural jacobian is given by: 


7 [9 ON ON 
oo Gu Os ales WO 


Device Jacobian 


The device jacobian is given by: 


ete —P.q  —Psp 
-“\0 0 =P, 0 
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D-2: Network Description 


Figure D-2-1 1s a block diagram of the system represented in Figure D-1: 


Figure D-2-1: Parallel Generator Example Block Diagram 





D-2.1: Network Specification 


Using the same variable labeling convention as in Appendix C, the devices and 


network are specified by: 


VDS Generator GI 


Terminal Potential Variable Flow Variable 
VQ XGI_V_nep XG1_Q nif 

DP XGI_D_nep XG1_P_nif 

p XG1_p_iep 

q XG1_q icp 

Parameters: 

Ve 1.05 PU 

Dg 0.00 rad 

P, 1.00 PU 

Import Vector: 
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Node 


NNW — 





PQS Generator G2 


Terminal Potential Variable Flow Variable Node 
vo x G2_V_nip XG) _O_nef | 
DP XG2_D_nip XG1_P_nef g 
Pp XG2 _p_tp > 
q XG2_q_iip 6 
Parameters: 
ge 0.50 PU 
Import Vector: 

XG2_V_nip V, 

Pa XG2_D_nip | __ V; 
G2_i = 
XG? p_iip Vs 
XG2_¢_iip Ve 


Transmission Line T12 


Terminal 
VQ!1 
DP1 


VQ2 
DP2 


Parameters: 
R 0 
xX (0. 


Import Vector: 


Potential Variable Flow Variable Node 
X712_V1_nip XT12_Q1_nef l 
X712_D1_nip X712_P1_nef 3 
X712_V2_nip X712_92_nef ju 
X712_D2_nip XT12_P2_nef 4 
PU 
PU 
XT12_V1_nip V, 
x | A712_D1_nip | _ V; 
Tina = 
XT}2_V2_nip V, 
x T12_D2_nip V, 
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PQ Load L3 


Terninal Potential Variable Flow Variable Node 
VQ X13 V_nip X13 O_nef 2 

DP X13 D_nip X13 _P_nef 4 
Parameters: 

lie 0.60 PU 

0, 0.10 PU 


Import Vector: 


ee X13 _Vnip | _ Vy 
cial XL3_D_nip Va 
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D-2.2: System Variables and Equations 


There are eight system variables and equations associated with this example. There 


are the six node potentials plus two import flow variables ordered in the following manner: 


T 
Xy= M1 Vo Vg Va Vas Ve lero lor P| 


The eight system equations are composed of four Kirchhoff Current Law equations 


and four potential equations: 
g 1(%,,5) = 16) 9 +%G2_¢9 nef +*712_01 nef 
82(X,,5) = X13 0 nef © *712_02 nef 
& 3(%5) =16, p+ XG2_p_nef *X712_P1_nef 
g 4(%,5) =Xp3 Pp nef +X712_P2 nef 
81_61_vXn5) = Vie XG1_V_nep 
83 G1_pXsys) = V3 —XG1_D_nep 
85 G1_pAsys) = Vs —XG1_p_iep 


g 6_G1_q OG) = Ve —~XG) _9_lep 
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D-2.3: System Structural Jacobian Matrix 


Using the device structural jacobian matrices along with the system equations, the 


following system structural Jacobian can be created: 


oo o™ 2 24224 24 
oooo 2222 
OF =e See 
oooo zzz 2 
SS S.o- oe eor 
+O 2 Close 2c 
ao 2 2 Ge as 
ee Oy 1S (ne 


Applying the system reduction algorithms, five blocks can be identified: two 1x1 


element blocks and three 2X2 element blocks: 
Block I 


System Row: 5 
System Column: 1 
System Variable: V;, 


Equation: 
81 G1_vX%ys) = Vi -Xe1_v rep 
Structural Jacobian: 


Jp, = (2) 
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Block 2 


System Row: 
System Column: 3 
System Variable: V; 


Equation: 
& 3. G1_D Xs) =V3—-XG1_D_nep 
Structural Jacobian: 
Jp2 = V7) 
Block 3 


System Rows: 2 4 


ho 
> 


System Columns: 
System Variables: V,  V, 


Equations: 
£,(x,,,) = X13 0 nef *XT12_Q2 nef 


8 4(Xsy5) = Xr3_p_neg + X112_P2_nef 


Structural Jacobian: 


N WN 
Ju=| y Hf 
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System Rows: 
System Columns: 
System Variables: 


Equations: 


Structural Jacobian: 


System Rows: 
System Columns: 
System Variables: 


Equations: 


Structural Jacobian: 


Block 4 


5 fi 
5 8 
Vs dep 


g 3(%,,,) = 161 p+XG2_p_nep +X112_P1_nef 


85 G1_pXays) a V, AGI ip iep 


| | 
J 34 =| 4 
Block 5 
] 
6 
V, I¢1 9 


84(% 5) = 1619 +XG2_9 nef + X12 Ql _nef 


§6Gi og pe) a Ve ~~ XG1 _q_iep 


Kaen 
Jns=|j a 
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D-2.4: Solving the System 


Applying the equations for the first two blocks yields: 


0 
V, = | } =1.05 PU 
0 


The remaining blocks are systems of 2X2 equations and unknowns. Blocks 3 and 5 
are nonlinear and must be solved iteratively. Block 4 is linear block requiring only one 


iteration: 


Block 3: 


aa 1.0000 -(). 0.0293 0.0085 
ows] sore] tases 





Block 5: 


| | we | w@ | tad | 


0. 0000| ).0000 0. 1750] 9.0000 


03828 aa00e.é] 000. 
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D-3: Summary of Results 


Bus 1 
Bus Voltage Magnitude 1 Oss 3) 
Bus Voltage Angle 0.00 rad 
G1 Real/Reactive Power -0.4125 PU -0.1167 PU 
G2 Real/Reactive Power -0.2063 PU -0.0583 PU 
T12 Real/Reactive Power 0.6188 PU 0.1750 PU 
Bus 2 
Bus Voltage Magnitude 0.9933 PU 
Bus Voltage Angle -0.1105 rad 
L3 Real/Reactive Power 0.6000 PU 0.1000 PU 
T12 Real/Reactive Power -0.6000 PU -0.1000 PU 


Information Node 5 (p) 
Magnitude 0.4125 
Information Node 6 (q) 


Magnitude 0.2828 


SES) = 





Appendix E: Waveform Examples 
E-1 Examples of Waveform Types 


While the possibilities of waveform definitions is endless, this thesis will concentrate 


on the following waveform types: 


Data Series 
Fourier Series 


Legendre Series 
Polynomials 


Matlab Polynomials 


Chebyshev Senes 





The code in the above table refers to the value of the type element in the WAVEFORM 


structure. 
E-1.1 Data Series 


A data series consists of m equally spaced samples of the waveform stored in an array 
of double precision floating point numbers. The first coefficient 1s associated with the 
value of the waveform at the beginning of the time interval and the last coefficient is 
associated with the value of the waveform at the end of the time interval. Each element of 


the array is given by: 


c;=f(t) 


j 
fe ye mae Ua 


A data series representation is primarily used for plotting the time history of 
variables and for calculating waveform operators which would prove difficult with other 


waveform types. 
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E-1.2 Polynomial Expansion 


A polynomial expansion consists of n coefficients of a polynomial representation of 


the waveform normalized over the interval [-1 1]. 





C= cn 
i=] 

[—-[ 
x=-1+2 

[,—[ 


Polynomial expansions are useful for evaluating switching operators described 


above. 
A Matlab Polynomial expansion is expressed in descending order: 


ie) = S ou 


i=] 


E-1.3 Orthogonal Function Series 


Orthogonal Function Series can be an excellent means for representing waveforms. 
In an orthogonal series representation, the value of the coefficient of a given order of the 
characteristic function 1s independent of the number of terms in the orthogonal series. This 
means truncating an orthogonal series by eliminating higher order coefficients will still 


result in the best possible fit with the remaining coefficients. 


In general, an orthogonal series representation is of the form: 


fQx)= x CF; _\() 
Ne 4, eek. | 
where Fx) is the ith order characteristic function of the orthogonal function series 
with respect to the weighting function r(x). These characteristic functions observe the 


following property: 
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Oo 





fre )\F_A(xX)F (x)dx =O for men 


‘a 


[ reoR QF (dy = G(m) 
Xo 


With this property, the coefficients c; of the series can be found: 


- 
‘” Gi-1) 





[ rearooF, .eoae 
*o 


E-1.3.1 Fourier Series 


Perhaps the most widely used orthogonal function series is the Fourier Series. 
Unfortunately, the Fourier Series is unsuitable for dynamic simulations. To see why, one 


need only look at the manner in which a function is expressed in a Founer Series: 


f(x) =A) + y A, cos(imx) + y B. sin(inx) 





t —fo 
x=-1+2 
f, oa lo 
Notice that at x=J and x=-J sin(inx)=0 and cos(inx)=(-1)' . Consequently 


f(1) = f(-1). In other words, the starting value and ending value of any waveform 
represented by a fourier series is forced to be identical. In dynamic simulations however, 


we often have equations of the form: 


a) 


= f(2: 
ay ay) 


This equation is normally evaluated by integration: 


y=Yt | te, oat 
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where y, is the value of the waveform y at the beginning of the interval. If y is 
represented by a Fourier Series, then y evaluated at the end of the interval will also be yy. 
In other words, while the value of a state variable may change within the interior of a 
time interval, at the boundaries, the value is constrained to be a constant independent of 
the length of the time interval. This constraint is artificial and not a property of real 
physical systems. 


E-1.3.2 Legendre Series 


Legendre Series use legendre polynomials to form the basis of an orthogonal 
function series over the interval [-1 1]. Legendre polynomials L(x) of order i are defined 


by the following equations: 


L u;(x) f ; 
A(x) = aa or 1 even 
coe ie 7 
(x)= va) or 1 Oo 


G41) 2 1E=DEF DEF) 4 = DE-MG+ DEE HS) 


ay 4! 6! 


(Siar 4) 3 n G@-1G-3)@+2)@ sie) eos 


a eT 5! 


ee ee 


The first six Legendre polynomials are readily found to be: 
L(x) = 1 


EO) 


] 
La) =5 Gx" - 1) 


L,x)= = (5x" — 3x) 
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1 
L(x) = 5 35x" — 30x° +3) 


L(x) = = (63x° ~70x? + 15x) 


Legendre Series also obey the following recursion formula 
GH Oo) = Cat xb) =n) 


An nth order legendre series representation of a waveform is given by: 


f(x)= x c,L;_ (x) 





where: 
X=-1 
x,=1 
r(x)=1 
F(x) =L() 
Seas 


The time interval [f,4,] can be mapped to the interval [-1 1] with the following 
transformation: 
t abe 
f, ae lo 


X= —l +2 





The coefficients c; can be found by integration: 
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(= J fot, a 


E-1.3.3 Chebyshev Series 


Chebyshev Series use Chebyshev polynomials to form the basis of an orthogonal 
function series over the interval [-1 1]. Chebyshev polynomials T(x) of order i are 


defined by the following equations: 
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T(x) = 1 
T(x) =x 
PO ee = eX )retor f 23 
The following three Chebyshev polynomials are given by: 
TC) es | 
T,(x) = 4x° -— 3x 
T, (x) = 8x* ~8x7 +1 
An nth order Chebyshev Series representation of a waveform is given by: 
f(x) a 2 CT; - (x) 
where 
k= a! 
Da 


the weighting function r(x) is given by: 





rx)=— 
1-x? 
and: 
F(x) = Tx) 
G(O)=n 


G(m) => for m>0Q 
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E-2 Waveform Conversions 


This section describes how to convert a waveform consisting of a vector of 
coefficients of order n, to a waveform of possibly a different type composed of a vector of 
coefficients of another order n,. In all cases, the conversion is a linear matrix operator. 


Hence for given values of n, and n,, the conversion matrix need only be calculated once. 


From here on, L(x) refers to a vector containing the polynomial coefficients of the ith 
order Legendre Polynomial. L{x;) refers to the ith order Legendre Polynomial evaluated at 
x, Likewise, T;(x) refers to a vector containing the polynomial coefficients of the ith order 


Chebyshev Polynomial. 7(x;) refers to the ith order Chebyshev Polynomial evaluated at x,. 
E-2.1 Legendre Series 
E-2.1.1 Legendre Series to Data Series 


Converting a Legendre Series of order n, to a data series of order n, requires the 


construction of the following matrix: 





1 L£,(Q%) Li(%) + EN ~1(%) 

ib GC) LCS ae eum ein Cay 

eee MeSCa a Los 

Arp = 

I Li, -1) LQ, -1) ae tee) 

l 
xX =—1 +2 

n,—1 


If C, is the vector of the Legendre Series coefficients and Cz is the vector of data 


senes points, the following relation holds: 
C,=A,,C, 
E-2.1.2 Legendre Series to Legendre Series 


Converting a Legendre Series of order 2, to order 7, requires only the truncation of 


terms if n, > , or the insertion of zeros in the higher order terms if 7, < nm. 
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E-2.1.3 Legendre Series to Chebyshev Series 


Converting a Legendre Series of order n, to a Chebyshev Series of order 7, first 
requires the truncation or padding with zeros of the Legendre senes to order n,. The 


resulting Legendre Series should then be multiplied by the following upper trianglular 


matrix: 
7 mee an 

Ae ame) EK) ci) Tee) 

pete) Tx) TL) Ro De) 


where L(x) is a vector of order n, holding the polynomial coefficients of the ith order 
Legendre Polynomial and Tx) is a vector of order n, holding the polynomial coefficients 


of the ith order Chebyshev Polynomial. 
E-2.1.4 Legendre Series to Polynomial Expansion 


Converting a Legendre Series of order n, to a polynomial expansion of order n, first 
requires the truncation or padding with zeros of the Legendre senes to order n,. The 
resulting Legendre Series vector should then be multiplied by the following upper 


triangular matrix: 
IN, = (ECoG a Ao Murr Ames (69) 


where L{x) is a vector of order n, holding the polynomial coefficients of the ith order 


Legendre Polynomial. 
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E-2.2 Chebyshev Series 
E-2.2.1 Chebyshev Series to Data Series 


Converting a Chebyshev Series of order n, to a data series of order nm, requires the 


construction of the following matrix: 





ie FO) T,%) «oT, 1%) 
ay ee Ty%) ee, 4) 
] Sees) 10) — i es) 
Arp = 
T(x, -1) T(%,,-1) =< jE ea ern 
= — le 
Ny a 


If C, is the vector of the Chebyshev Series coefficients and C, 1s the vector of data 


series points, the following relation holds: 
Cy = Arp C, 
E-2.2.2 Chebyshev Series to Legendre Series 


Converting a Chebyshev Series of order n, to a Legendre Series of order n, first 
requires the truncation or padding with zeros of the Chebyshev series to order n,. The 


resulting Chebyshev Series should then be multiplied by the following upper trianglular 


matrix: 
A= A Ae 
A, =o) LiG@) Lx)... LAG] 
Pea) Gy BG) a5 120) 


where L{x) is a vector of order nm, holding the polynomial coefficients of the ith order 
Legendre Polynomial and T{x) is a vector of order n, holding the polynomial coefficients 


of the ith order Chebyshev Polynomial. 
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E-2.2.3 Chebyshev Series to Chebyshev Series 


Converting a Chebyshev Series of order n, to order n, requires only the truncation 


of terms if m, > 1, or the insertion of zeros in the higher order terms if 1, < 7. 
E-2.2.4 Chebyshev Series to Polynomial Expansion 


Converting a Chebyshev Series of order 2, to a polynomial expansion of order n, 
first requires the truncation or padding with zeros of the Chebyshev series to order 7. 
The resulting Chebyshev Series vector should then be multiplied by the following upper 


triangular matrix: 
A; = (ES) HC 56) eae Hes i(X)) 


where T(x) is a vector of order n, holding the polynomial coefficients of the ith order 


Chebyshev Polynomial. 
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E-2.3 Polynomial Expansion 
E-2.3.1 Polynomial Expansion to Data Series 


Converting a polynomial expansion of order n, to a data senes of order n, requires 


the construction of the following matrix: 





> n,—] 
1 Xo Xo ‘i 
9. n,-1 
ex, Xe 5 te 
. ny-l 
I) xy x4 x4 
App = 
2 n=l 
ea eee 
; 
Ny — 


If C, is the vector of the polynomial coefficients and C, is the vector of data senes 


points, the following relation holds: 
Cy =AppC, 
E-2.3.2 Polynomial Expansion to Legendre Series 


Converting a polynomial expansion of order n, to a Legendre Series of order n, 
requires first converting to a Legendre series of order n, then converting the Legendre 
Series to order m,. Recall that the matrix for converting from a Legendre Series to a 
Polynomial is upper triangular. Hence one only needs to use backward substitution to 


solve for the Legendre Series coefficients: 
A, = [Lo(*) Lx) £4)... L,-1X)] 
C,=A,C, 


E-2.3.3 Polynomial Expansion to Chebyshev Series 


Converting a polynomial expansion of order n, to a Chebyshev Series of order n, 


requires first converting to a Chebyshev series of order n, then converting the Chebyshev 
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Series to order m,. Recall that the matrix for converting from a Chebyshev Series to a 
Polynomial is upper triangular. Hence one only needs to use backward substitution to 


solve for the Chebyshev Series coefficients: 


A, = (To) ATO Xe 3 I, -14)) 


C,=A,C 


P 


f 


E-2.3.4 Polynomial Expansion to Polynomial Expansion 


Converting a polynomial expansion to another polynomial expansion of higher 
order only requires setting the higher order terms to zero. Converting to a lower number 
of terms requires more effort. The best method is to convert to an orthogonal function 
series, truncate, and convert back. Since all of these operations are linear matrix 
operations, the conversion matrix need only be calculated once. For this conversion, 
either the Legendre Series or the Chebyshev series would be appropnate since the type 
conversions to and from the series solution does not add any truncation error (The 
truncation error is solely due to the truncation of the Legendre Series or Chebyshev senes 


and not due to the conversions). 
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E-2.4 Data Series 


E-2.4.1 Data Series to Data Series 


There are many methods for converting a data series to another data series with a 
different number of coefficients. Two common interpolation schemes for performing this 
conversion are linear interpolation and cubic splines. These methods can be found in 


many numerical methods textbooks and will not be described here. 
E-2.4.2 Data Series to Legendre Series 


If n,2n,, a Data Series can be converted to a Legendre Series by taking the 
pseudo-inverse of the matrix converting a Legendre Series to a Data Series. If, <n), the 
Data Series can be converted in a similar manner to a Legendre Series of order n, padded 


with zeros to order np. 





i eG) Lj%) ss L -(%) 
t° 2£,6,) L(x) eB 4y) 
Po L,0Q) £0) ee bg 15) 
Arp = 
| L(x, -1) L(x, -1) vee fey, (Xp. 1) 
oe led 2 : 


n,-1 


If Cz is the vector of data series points and C, is the vector of Legendre Series 


Coefficients, the following relation holds: 
Cy = ApS, 
Ci=(AipAw) AipCy 
EK -2.4.3 Data Series to Chebyshev Series 


Converting a Data Series to a Chebyshev Series can be done in the same manner as 


the conversion to a Legendre Series: 
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| LOG) T (Xo) nee if _(%) 





] T,(x,) Ce) ee i =a) 
LO Qs) Ty(%,) «Ti -,G) 
Arp = 
1 T(x, -1) T(x, -1) vee T,, -1%n, -1) 
1 
x,=-1+2 
ny aa 


If C, is the vector of data series points and C, is the vector of Chebyshev Series 


Coefficients, the following relation holds: 


C,= ArpC, 
-1 
C.= (ArpAzp) ees 


E-2.4.4 Data Series to Polynomial Expansion 


Converting a data series to a polynomial expansion of equal or less order using a 
least squares fit is a straight forward process. If the number of points in the data series n, 
is equal to or Jess than the number of points in the polynomial n,, the resulting 
polynomial will pass through each point of the data series. If larger, the polynomial will 
not necessarily pass through all of the data series points, but will be a least square 


approximation. 
Forn, <n,: 
For this case, the problem is to solve for the coefficients of the polynomial c,; for 


i<n,. For the higher coefficients (> ,), c,;=90. In the following discussion, let C, be 


the vector of nm, data series coefficients and C, be the vector containing the first n, 


polynomial coefficients. Define the n, < n, matrix A as follows: 
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fe Ae 

a3 Be ae 
2 ani 

] x4 Ay 2 

App = 
> ny-l 
] Maen ne ce 
1 


ny—l 


Matrix App 1s Square, clearly has rank n,, and therefore 1s invertible. Consequently 
solving for C, is straight forward: 
Appt, = Cy 
C, = Ane, 


Forn,>n,: 


If the number of data points is greater than the number of polynomials, the number 
of columns in the App) matrix described above would have n, columns and n, rows. App 
would clearly not be invertible. The pseudo-inverse of App can be calculated and 


provides the least squares fit of the data series: 


-1 
Ce (AppApp) ApnC, 
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E-3 Waveform Arithmetic 


This section descnbes how to perform addition, subtraction, multiplication, and 


division on the various types of waveforms. 
E-3.1 Data Series 


Performing waveform arithmetic on data senes is very easy. The waveforms are 
converted to the proper size and then added, subtracted, multiphed or divided element by 


element. 


E-3.2 Polynomials 
E-3.2.1 Addition/Subtraction 


Adding or subtracting two polynomial waveforms simply entails converting the two 


waveforms to the proper length and adding or subtracting element by element. 
E-3.2.2 Multiplication 


Multiplying polynomial waveform W of size n, and Y of size n, together to get 
polynomial Z of size n,+n,-1 can be accomplished by constructing the following 


matrix of sizen, +n,-1Xxn,: 


ar OO 0 0 

Y, Y¥, 0 O 0 0 

Y, Y, Y, 0 0 0 

eee jas 
M,= 

0 0 Yeni 

0 0 Y¥ 


Z=M,W 


Z can now be truncated or padded with zeros to convert it to the proper length. The 


truncation of a polynomial is discussed in section E-2.3.4. 


Oa 





E-3.2.3 Division 


Dividing two polynomial expansions can be difficult, particularly if the 
denominator polynomial has one or more zero crossings. In general, there 1s no simple 
method for performing the division, although the recursion process described in this 


section will work. Define the problem to be: 


S Bxi7? 
a re ie 


I 
1 me 


ee 


There are two parts to the problem. The first task is to use synthetic division until 
the numerator of the remainder is of size n,-1 or less. The second task is to convert the 
remaining fraction into another polynomial expansion by a process similar to synthetic 


division, but proceeding from the constant term and working up in order. 


Synthetic division is the process of dividing one polynomial by another until the 


remainder is of order | less than the denominator. 


ae l=] [-1 
Dak a »y VX 
a es oe? as =] 
4] ia n 
c Cc C 
k=] Ly k-] 
Dae pare 
k=] k=] 
Ma 
r,=d,-—c¢, 
ne. 
d,, 


¥fing-n, +1) a 6. 


Initially, d; is set equal to b;. After the first iteration, d; is set equal to the remainder 
r, The process is repeated until nz =n,-1. At this point, the direction of the division is 


reversed and we get: 


ny n= 1 

fad fA 
» dx d Dar x 
[= 1 |=} 
n ae. n 
¢ Gc c 

k-1 1 k-1 
DEG Dex 
k=] k=} 
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ee 
Pop d= = C 
{—] l Cy l 

In this manner, we can express the remaining fraction as another polynomial 
expansion. The actual values for Y; are equal to the sum of the components from the 


forward and backwards synthetic division. 


Note that if the denominator has a zero over the interval [-1,1], the backwards 


synthetic division will result in a diverging series. 
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E-3.3 Legendre Series 
E-3.3.1 Addition/Subtraction 


Adding or subtracting two Legendre Series waveforms simply entails converting the 


two waveforms to the proper length and adding or subtracting element by element. 
E-3.3.2 Multiplication 


Multiplying two Legendre Series together can be accomplished in two ways. The 
first way is to convert the Legendre Series to polynomial expansions, multiply the two 
together, then convert the product to the Legendre Series of the proper size. The second 


method uses the recursion formula for the Legendre series to assist in the process: 


To multiply Legendre Series Y of size n, by the Legendre Series W of size n,, to 
obtain the Legendre Series Z of size n,=n,+n,,-1, Y must first be converted to a 


polynomial expansion Y, of size n,: 
A, =[Lox) £4) £0) ... L,, 1)] 


Y,=A,Y 


P 


The recursion formula for the Legendre Series is given by: 


4] 
xL(x)= [J e+( i Jeno 


which can be translated into the following matrix for multiplying a given Legendre 





Series of size n, by x: 
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Teas 0 0 
3 
2 
eee 0 0 
5 
2 3 
nee eeon 0 0 
3 7 
3 
OuEO= =O 0 0 
5 
ee 
00 0 0 cia 0 
POG oy 
000 0 0 oe 
i Rees 
00 0 0 ie 0 
= =e 


If we define the vector Y p1 tO beny. P padded with zeros such that it is of size n,, we 


can define the following n, X n,, matrix: 


Pee a Vea Ae eay | AREY 


p! p! p! p! 
The final n, x n, multiplication matrix A,,, can now be found: 


Ann = AmpiAr 


mpl 
Za AW 


Of course Z may have to be truncated or padded with zeros to convert it to the 
desired length. 
E-3.3.3 Division 

There is no straight forward method for dividing two legendre series. The easiest 


way appears to be converting to polynomial expansions, performing the divison, then 


converting back to the legendre series. 


2454 


E-3.4 Chebyshev Series 
E-3.4.1 Addition/Subtraction 


Adding or subtracting two Chebyshev Series waveforms simply entails converting 


the two waveforms to the proper length and adding or subtracting element by element. 
E-3.4.2 Multiplication 


Multiplying two Chebyshev Series together can be accomplished in two ways. The 
first way is to convert the Chebyshev Series to polynomial expansions, multiply the two 
together, then convert the product to the Chebyshev Series of the proper size. The second 
and preferred method uses an alternate definition of a Chebyshev Polynomial to assist in 


the process: 
T,,(x) = cos(n cos”'(x)) 
TAO LC) 


From this definition, the product of two Chebyshev Polynomials can easily be 


derived: 


T(x)T_, (x) = cos(n cos '(x)) cos(m cos '(x)) 
T,(x)T__(x) = : (cos((n +m) cos '(x)) +.cos((n —m) cos '(x))) 
TIT (2) =5 Tx on (0) +E yn) 
To multiply Chebyshev Series Y of size n, by the Chebyshev Series W of size n,, to 


obtain the Chebyshev Series Z of size n,=n,+n,,- 1, three n, Xn, matrices should first 


be constnicted: 
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NG 3 
_ 
Ne 
y, 0 0 0 0 
Y, Y, 0 0 0 
Y, ¥, Y, 0 0 
Aan a 
0 Y,,, ee 
0 0 ee 
Y, Y, Y, 
Y, Y; Y, 
ae - "3 "4 3 
0 0 0 QO 
Uae fa Poe 2 
0-05 2 OY, 
Jane =|0 O 0 ie 


The final m, X n, multiplication matrix A,,,, 1S given by: 


mi2 


l 
ee == 9 (Ars a; A 7 Ain3) 


Z=A,_,W 


- 247 - 








Of course Z may have to be truncated or padded with zeros to convert it to the 


desired length. 
E-3.4.3 Division 


There is no straight forward method for dividing two chebyshev senes. The easiest 
way appears to be converting to polynomial expansions, performing the divison, then 


converting back to the chebyshev series. 
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E-4 Waveform Functions 
E-4.1 Data Series 
E-4.1.1 Trigonometric and Exponential Functions 


All trigonometric and exponential functions can be performed point by point on the 


data series coefficients. 
E-4.1.2 Integration and Differentiation 


There are a number of techniques for integrating or differentiating Data Series. All 
are by their nature approximations and can suffer from numerical instability problems 


associated with conventional simulations. One simple method of integration employs the 





trapezoidal rule: 
U0 OR 6 ne OF? 2020 
h oh 
= = 0 0 O 
ee 
h h 
— hk, — 
7 Oo = SOF 
h h 
i oe ae 
5 5 8 Or O30 
a at, ye 0 O O 
yi 2 
Spp 
h h 
= oa 
5 h ie 0. 
h h 
5 tf Ah Jie ee 
h h 
ee ee (elias 
pee 
n-1 


With this matrix, the integral equation: 


- 249 - 





) i» eg { Wdt 


te] 
Becomes the matrix operation: 
Y=S,)>W+Y, 
The vector Y may be now converted to a different length if so desired. 


Differentiating a Data series can be done in a number of ways. The secant method 


can be easily implemented with the following matrix: 





—] 1 0 OO 0. ... 0 0 O 
1 ] 
—=— — es 0 
5 0 5 0 0 0 O 
1 1 
0 ar 0 5 0 0 0 O 
] 1 
1} O 0 —- 0-= 0 0 O 
SO) 2 2 
] ] 
0 0 0 0 0 -5 0 5 
0 0 0 O 0 0 -1 1 
] 
aaa 


Another approach 1s to choose a differentiation matrix such that it 1s consistent with 


the integration matrix. Consistency is defined by the following matrix equation: 
SppApp2 =M 


where 
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20 OSU 
Sal OO 
=i 0 1-0 
Veter Vie" 


The M matrix reflects the fact that differentiating a data series will destroy the 
subsequent constant of integration. Since Spy 1s generally singular, only its 


pseudo-inverse can be taken: 
T lar 
Dyp2 = (SppSpp) SppM 
This matrix actually has a very simple construction: 


2n —3 2n-5 2n-7 2n —9 
































d 
: n n n n 
] 2n —5 2n-7 2n-9 
nh ae n ty, n 
l 3 2n—-7 2n —9 
n “n ds; n ee 
1 
ioy ae ae 3 5 2n—9 
S| = == di, 
n n n n 
3 ] 
ae S — d 
n n n n 2 


d 


XX 


is equal to the the negative sum of all the other terms in the xth row. As a 


consequence, all of the row sums of Dpp, are equal to zero and Dpyp, is singular. 
E-4.1.3 Switching Functions 


All switching functions can be performed point by point on the data series 


coefficients. 
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E-4.1.4 Waveform Smoothing 


There are times when it may be desirable to remove high spectral content features 
of a waveform. One way to do this is to replace the value at each point in the time 
domain by the average of the waveform over some interval [x - A,x + A]. This can be 


accomplished by defining the following: 


nf =D) 
n, = in 
2 


where int(x) is the integer nearest x. 






































— = = 0 0 

Wn ne 7h Dn 

1 1 ] l 1 0 
nyt+1 ntl nytl n+1 nytl 

1 1 1 1 1 1 
Net2 Mmet2 mgt2  mgt2 mt2 nyt+2 

Ao 
1 1 1 1 1 1 
. ee a ee 
Die 5) ne on, a. 


Multiplying a data series by A,,,,, Will return a smoothed version of the data series. 
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E-4.2 Polynomial Expansion 


E-4.2.1 Trigonometric and Exponential Functions 


There is often no direct way of evaluating a trigonometric or expontential function 


of a polynomial expansion. Instead, the function is performed on a data series converted 


from the argument polynomial. The resulting polynomial is then reconverted back into a 


polynomial. 


E-4.2.2 Integration and Differentiation 


Integrating a polynomial Y of size n, results in another polynomial Z of size 


n,=n,+1. Then, xn, integration matrix Spp is given by: 


eds 

2 

1 0O 

] 

9 3 

0 O 

Spp=| 0 0 
0 O 

0 O 


The integral is evaluated by: 





ee. 

3° 4 noe 
0 O 0 

0 O 0 

] 

- QO 0 

3 

1 
0 4 0 
0 O 
no) 

OO 0) 
L= Sp) +Z, 





Of course, Z may be converted to a polynomial of a different size if desired. 


Differentiating a polynomial Y of size n, results in another polynomial Z of size 


n,=n,-1. The n,x n, differentiation matrix App 1s given by: 
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0 1 0 0 0 0 
Oe. 2 0 0 
OUI Os 0 0 
App 
0 ia 0 
0 0 n,-1 


The Differential is evaluated by: 


Z =AppY 
E-4.2.3 Switching Functions 


Switching functions are those which produce a Polynomial waveform Y which 1s 
composed of m pieces of other Polynomial waveforms. Let f; be the Polynmomial 
representation of the jth piece of Y. Let x,(y) be the x coordinate of the ending point of 
the jth piece where x,(0) = -1 and x,(m) = 1. 


Define Y, to be the Legendre Series representation of Y : 


Then using the orthogonality property of the Legendre Series: 


Xo) 


{ 2 - Jjotaener 


Xp - 1) 


Y,, = 


‘ 





We 


= 


J 


Now define the following row vector: 


I(x) = Lo) Li%G)) L200)» Ly -1QoG))  £,009))] 
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With Sp, as defined in section E4.3.2 and A,,,4) as defined in section E-3.3.2 the 


solution for Y, can easily be found: 


1 
- 0 QO 0 
Z 
3 
05 8 0 
5 
0 0 = 0 
G,= ye 
2n—\ 
000. Met 


m 


ee x (I) - HU — I) SpAmpi(f)G, 


=] 
Now we need only convert Y, to a polynomial exansion Y: 


Y=A,Y, 


where 
A, = (Lox) L(x) L(x)... L,, -(*)] 


E-4.2.4 Waveform Smoothing 


There are times when it may be desirable to remove high spectral content features 
of a waveform. One way to do this is to replace the value at each point in the time 
domain by the average of the waveform over some interval [x - A,x + A]. This can be 


expressed by the following integral: 


x+A 
n ; l n y 
Y.x° = Wr o'dt 
2 mn 2A Py y 
x-A 
The only problem with the above equation is near the boundaries x = -1 and x = I 


where the integration interval has the possibility of crossing the boundaries and including 
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within the average a section of the polynomial outside the defining interval [-1,1]. Hence 
the smoothed polynomial should be composed of the following three segments (assuming 
= 1): 


~-l<x<-l+A 


x+A 


n 1 
oo Wri 'dt 
2 ai x+At+1 J 2 
—-l+Asx<l1-A 
1 x+A 
Yx'7 == | Wt dt 
py De & : 
x-A 
1-A<x<l 
SY f3 Wee ‘dt 
2 a a i 2 


Note, if 1 < A <2 then the interval boundaries are given by: 


[(l1-A , -1+A] 
[(-l+A , 1] 


If A> 2 then there is only one interval and the average of the waveform is retumed: 
1 1 
=5{ > W,x' "dx 
2 aff i=] 


Y=) On we 


For A < 2 evaluating the integrals require the definition of shifting a waveform left 


or right by A. This can be done by constructing the following binonomial matrix: 
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eee So 
a20O0- 
2OfF Ne 
SS a 
KH LAP 


exp 


This matrix can be generated by the following recursion formula: 
Bs A1,J) = 1 
B.,(2:n, 1) =0 
By) oy tee =) 1) 


Multiplying B,,, element by element by the following matrix B, will give us the 


transformation matrix B,,,, for shifting a waveform left by A. 


owe A A 

OF eA eA ee: 

00 1AfnM 
ee 00 0 1 =A 

00 0 0 1 
x Ma = LN + Ayo! 
= = 


M =BayN 


The tools are now all present. W can be integrated using the integration matrix Spp 
defined in section E-4.2.2. The limits of integration for the three segments can be applied 
by either using B,,,, for the limits involving x, or by direct evaluation for those limits not 


involving x. Dividing by the averaging interval comes next. For the first and third 
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intervals, the methods outlined in section E-3.2.3 can be used to divide a polynomial by 
another polynomial. Finally, the procedure for generating Switching Functions described 


in section E-4.2.3 can be used to generate the coefficients for the solution Y. 
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E-4.3 Legendre Series 
E-4.3.1 Trigonometric and Exponential Functions 


There is often no direct way of evaluating a trigonometric or expontential function 
of a Legendre Series. Instead, the function 1s performed on a data series converted from 
the argument Legendre Series. The resulting polynomial is then reconverted back into a 


Legendre Series. 
E-4.3.2 Integration and Differentiation 


Differentiating a Legendre Series can be done easily by differentiating the recursion 


formula for the Legendre Series. Recall: 


(i+ IL, (4) = (21 + Dxh(x) -iLy_(%) 














Differentiating: 
Geog Qi ft wai) i \dL,_,(x) 
er a) | ee | sss FO.) (reed Pecan af Neceaeerrrieics 
ax [ttl dx i+] ax 
where: 
dL (x) 
ea 
dL,(x) 
a 1=L,(x) 
The goal is to generate the following n X n matrix: 
ae AL(x) dbx) dL,(x) aL, _ (x) 
eee ax dx dx Fes dx 


The columns of Ap, can be solved recursively once we define the matrix A,, for 


multiplying a Legendre Series Vector by x. 


: ne 
cL) = ( sic fen +| ft, 0) 
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ee ee 0 0 
2 
10 = 0 0 0 
2 
0505 0 0 
3 
00 = 0 0 0 
a= 
=) 
000 0 Tree ee 
4 
(urOmtOs (0:reae 0 an 
=| 
00 0 0 Taree 


Note that the last row of Ay, has been elminated to make the matrix square. This 
will not cause any problems since in the recursion formula, the last coefficient of the 


vector multiplying Ay, 1s always zero. 
Let Ap, (2) represent the jth column of Ap,. Let I be the n x n identity matrix. The 
recursion formula states: 


i 


2i+1 
au JAaAouts toh ath Geet) )) {outa 


i+] 





Ap (:,i + 2)= ( 


l<isn-2 


Once Ap, is constructed, it can be used to calculate derivatives. Let W and Y be 


vectors of Legendre Series coefficients of size n. Then the following statements are 


identical: 
dWw 
Y =— 
ax 
Y=A,,W 


Of course, the nth coefficient of Y will always be zero since the nth row (as well as 


the first column) of Ap, wul always be populated with zeros. 
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Integration 1s a bit more complex. In general, the problem 1s to solve the following 


equation: 
=a | Wat 
t=—] 


First, the 2+1 x n indefinite integral matrix S,, should be found. The easiest way of 
generating S,, begins by adding an additional column to Ap, using the same recursion 


formula to form the n x n+1 matrix Ap,,. S,, is simply the pseudo-inverse of Ap,;: 


-1 7 
Sip = GemAci) Apu 


The next step is to evaluate the integral at x =-1. This can be done by multplying 


the following row vector by S,,: 
Ka (lee tee ee sd (1) 
S_, = X_ Siz 


The first row of S,, contains all zeros. If this row is replaced by -S_, and the 
resulting matrix called S,,, we have all the pieces for calculating the integral of a 


Legendre Series: 
Y=5,, 1 +1, 
Of course, the vector Y may have to be truncated or padded with zeros as required. 


E-4.3.3 Switching Functions 


Switching functions are those which produce a Legendre Series waveform Y which 
is composed of m pieces of other Legendre Series waveforms. Let f; be the legendre 
series representation of the jth piece of Y. Let x,(y) be the x coordinate of the ending 


point of the yth piece where x,(0) = -1 and x,(m) = 1. 


Let: 
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Then using the orthogonality property of the Legendre Series: 


Xo) 


Y= > { = |p 
1°" 9G -)) 


Now define the following row vector: 


Fon) — 14 Cce)). e 60D) Gn) <4 O90) = GU) 





With S,, as defined in section E-4.3.2 and A,,4) as defined in section E-3.3.2 the 


solution for Y can easily be found: 


l 
- Q0 OO. ... 0 
Z 
3 
0 5 0 0 
5 
0 0 = 0 
G, = 2 
2n—1 
0 0 O ———_— 
Ds 


¥7 = ¥ UD) — HF W)SoAmul 6G; 


E-4.3.4 Waveform Smoothing 


There is no obvious method for performing waveform smoothing in the Legendre 
Series spectral domain. Instead, the waveform should be converted to a polynomial 


expansion and the methods of section E-4.2.4 employed. 
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E-4.4 Chebyshev Series 
E-4.4.1 Trigonometric and Exponential Functions 


There is often no direct way of evaluating a trigonometric or expontential function 
of a Chebyshev Series. Instead, the function is performed on a data series converted from 
the argument Chebyshev Series. The resulting polynomial 1s then reconverted back into a 


Chebyshev Series. 
E-4.4.2 Integration and Differentiation 


Differentiating a Chebyshev Series can be done easily by differentiating the 


recursion formula for the Chebyshev Polynomials. Recall: 


Ty 4(x) = 2xT (x) - T;_ (00) 





Differentiating: 
aT;,,(%) aT (x) aT; _ (x) 
ee a) 
dx Fe ee eare 
where 
dT(x) 
= 
dT,(x) 
as 
The goal 1s to generate the following n X n matrix: 
A =| tot) ati) T(x) aT, _\(X) 
AEE: dx ae, Bleeds 


The columns of Ap; can be solved recursively once we define matrix A,; for 


multiplying a Chebyshev Series vector by x. 


] 
Ta) = F(T) +) 
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It 
05, 9 0 eed) 
eet 0 O 
2 
l ] 
05 9 5 0 0 
00 = 0 0 0 
Z 
Ay; = 
: 0 
U0 O 6 5 
1 
ES OS 9 ae 8 0 5 
oie 
OOo ="0 5 


Note that the last row of Ay; has been eliminated to make the matrix square. This 
will not cause any problems since in the recursion formula which follows, the last 


coefficient of the vector multiplying A,; is always zero. 


Let Ap7(:y) represent the jthe column of Ap; Let I be the n Xn identity matrix. 


The recursion formula states: 
Apr.t + 2) = ZA,-A,,,(.51 + 1) 21,7 +1) —A,-C,1) 
1<isn-2 


Once Ap; has been constructed, it can be used to calculate derivatives. Let W and Y 


be vectors of Chebyshev Series coefficient of size n. Then the following statements are 


identical: 
dw 
Y =— 
dx 
Y =A,,;W 


Of course the nth coefficient of Y will always be zero since the mth row (as well as 


the first column) of Ap; will always be populated with zeros. 


yeu 





Integration is a bit more complex. In general, the problem is to solve the following 


equation: 


First, the n+1 x n indefinite integral matrix $,; should be found. The easiest way of 
generating S,,; begins by adding an additional column to Ap; using the same recursion 


formula to form the n X n+I matrix Ap;;. S;7 1S simply the pseudo-inverse of Ap7;: 
rg Sar 
Sr =(AprApri) Apri 


The next step is to evaluate the integral at x =-1. This can be done by multplying 


the following row vector by S,,: 


Gp, = (10) et ee (9 ie! 


See eng = 


The first row of S,,; contains all zeros. If this row is replaced by -S_, and the 
resulting matrix called Sp,, we have all the pieces for calculating the integral of a 


Chebyshev Series: 
Sn ie 
Of course, the vector Y may have to be truncated or padded with zeros as required. 
E-4.4.3 Switching Functions 


Switching functions for Chebyshev Series can not be evaluated as easily as the 
switching functions for Legendre Series due to the weighting function r(x) for the 
Chebyshev Polynomials. Recall: 


1 
] a 

C, =-{ fe) _ =aXx 
Jed Vi-x 


1 
2 (ee 
¢, == | ax 
es 


Te 


The situation is not hopeless due to the following integral equations: 
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| hee) 


jaye 





{ as 7 NS 


iene 








2m ae 
x (2m)! ee FU Geo)! ceeat mest \s) 
dx = VNl-x? Yo se 4 
Iz —? (m | pale 2 Oe oe 


2m +1 = 0) ie 
{4 PE eR aU 


1—x? r=1(2m +1)rty 





Thus if f; is the Chebyshev Series representation of the yth piece of Chebyshev 
Series waveform Y and x,{(j) is the x coordinate of the ending point of the jth piece, then 


we can State the following: 


Y, 
Y, 
Y, 
Vas 
Ge 
x,(0) =-1 
Xo(m) = 1 
“ef (x) 
| ps AX 
a= = > dy 
el as ee 


Xo) 
y a2% f ACTn@ 


" Tjat aN 
Xo — 
The process should now be clear: 


1. Convert f(x) to a polynomial representation f,{x) 
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2. Multiply f,{x) by the polynomial representation for 7;(x) and call the resulting 


PJ 


polynomial f,(x). 


3. Use the above integral equations to evaluate at x =x,(j) and x =x,(j-1) the 


integral of f,{x) term by term to form the yth component of Y; called Y;,. 
4. Sum up Y;, over j to produce Y;. 


While the above process will produce the correct values for Y,, the following 


method is much easier to calculate and produces nearly identical results: 
1. Convert f(x) to a Legendre Series representation f,{x) 


2. Calculate the Legendre Series Representation Y, of Y with the methods of section 
E-4.3.3. 


3. Convert Y, to the Chebyshev Series Representation Y. 


E-4.4.4 Waveform Smoothing 


There is no obvious method for performing waveform smoothing in the Chebyshev 
Series spectral domain. Instead, the waveform should be converted to a polynomial 


expansion and the methods of section E-4.2.4 employed. 


- 267 - 








Appendix F: Model Development 


The following electrical power system models have been develped in support of 
WAVESIM: 


Three Phase Synchronous Generator 
Voltage Regulator 

Prime Mover 

Three Phase Switch 

Transmission Line 

Constant Impedance Loads 
Reduction Gear 

Propeller 

Ship Dynamics 

Pulse Generator 


Induction Motor 
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F-1 3 Phase Synchronous Machine Model 


Two models are presented for simulating a three phase synchronous model. The first 
expresses the voltages and currents in terms of a rotating reference frame (dqQ) rotating at 
the base frequency. This model is suitable for studies where the voltages and currents are 
balanced, nearly sinusoidal, and near the base frequency. For fast transients or unbalanced 
operations, the actual instantaneous values for the voltages and currents should be used (abc 
frame). Both models are very similar in that the terminal values are transformed to a 


rotating reference frame alligned with the rotor of the machine (Park’s Transformation) 


F-1.1 DQO Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Stator Direct Vp (import) Ip (export) (1) Normal 

Stator Quadrature Vo (import) Ig (export) (1) Normal 

Stator Zero Sequence V, Gmport) I, (export) (1) Normal 
Mechanical (,, (amport) T., (export) (0) Normal 

Field Voltage Vep (mport) Information 

Stator D-axis Current Ip (export) Information 

Stator Q-axis Current Tq, (export) Information 

Stator 0-axis Current I (export) Information 

Field Current Tq (export) Information 


The import x;,,, and export x,,, vectors are defined by: 


Vp ly 

Vo Io 

x imp - Vo I 
Ve D _ I F 

O,, | ae 

Lp, 

Lor 

ly 
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Parameters 


Xd 
Xq 
? 


X4 


ce 


Xg 


States 


Os 
Was 
Woes 


” 
€ gS 
C45” 
, 
é gS 


Synchronous Reactance (PU) 
Negative Sequence Reactance (PU) 
Transient Reactance (PU) 


D-axis Subtransient Reactance (PU) 

Q-axis Subtransient Reactance (PU) 

Annature Leakage Reactance (PU) 

Transient Open Circuit Time Constant (seconds) 
D-axis Subtransient OC Time Constant (seconds) 
Q-axis Subtransient OC Time Constant (seconds) 
Armature Time Constant (sec) 


Inertia Constant (sec) 
Pole Pairs 


Field Current for no load rated voltage (amps) 


Base System Frequency (rad/sec) 
Base System Angle (radians) 


Base System Voltage (volts) 
Base System Power (watts) 


Base Machine Voltage (volts) 
Base Machine Power (watts) 


rotor angle wrt to synchronous frame (rad) 

D-axis flux-linkage (PU) 

Q-axis flux-linkage (PU) 

Q-axis voltage behind subtransient reactance (PU) 
D-axis voltage behind subtransient reactance (PU) 
Q-axis voltage behind transient reactance (PU) 
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Equations 


Constant Definitions 


Base Quantities 





ee 
= 35 
_ se 

SB Osp 
I _ 2 me 
ro BV gs 
P 
[me MB 
Ops 


lg = pul Xa =a) 





Pup 
Veg =—— 
I 
pB 
Other Constants 
Xad ~ Xa %ar 
2 
oe 5 ae 
f Xa —XxXy 
Xf 
oes , 
(O,, 1 40 
a 
Xe = 
Xq ae ee 
Oe —x," 
a= , an 
Ag —~Xg 
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Angle Calculations 


O= fo. — 0,,P,)dt + Os 


0 =S(a,, - mPp) + Oso 


Ce = cos(O) 
Se = sin(9) 
Variable Rotation and Scaling 
Vy 
y 
v=] 7 
Vo 
Vid 
Vp 
V 
v=| ° 
Vo 
Vep 
V sp Vsp 
——M(C,.) -——Mi(S,) 90 0 
Vig S) Vie ( eC 
V V 
7 Me) G-M(Ce) 0 0 
MB MB 
K,= 7 
0 0 uses 
V up 
V 
0 0 ae 
Vp 
v=RV 
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Solving the electrical dynamical equations 


The five electrical dynamical equations must be solved simultaneously. Since the 
Integration Matrix and Multiplication Matrix are linear matrices, the entire problem 
becomes a linear process. Hence the system of equations can be represented by a matrix 


equation. 


First define the integration and multiplication matrices 
[xeoar =r 2 Set exe 


ee VEO eens Ve Sy 


Now we define the system of equations 











S S 
[+— —SM(@ —_ — 0 0 
T., (0,,P,) T., 
S S 
SM(@ [+— 0 —— 0 
(©,,P,) T., T., 
Be ” i ee S 
=| -—S§ 0 I+S 0 eee 
A i. if x , jap , Coe ay a? 
x,-x,” x 
0 S eee 0 [+S af 0 
Le Xq Uae Xq 
a—l 04 
0 0 —§ 0 [+S — 
fn T45 
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i) SO 0 420 
USO 3030 
0 0 7 0 0 
0007 0 
000 0 ¥F 


Bs 


| e325 
ea 3 
Soe a 
SSS ee 


So = 


BV eS, 


b= 


x=A'b 
= 4 





Calculating Export Variables 


First the currents in machine reference frame 


ly 
l 
ee 
lo 
Lig 
1 1 
ere 7 0 0 
Xg Xg 
‘= 0 as 
C= X q x q 
0 0 0 0 0 
x x 
0 0 = z 


ee 
Xaa(X¢ — Xta) Xaa(Xp— Xta) 


Now the currents in system reference frame 


ly 
I 
[= 2 
I 
lep 
I I 
—“M(C.) —M(S.) 0 0 
Isp Isp 
! Mp I up 
eS) SMG) 0 
R, Isp is Isp ? 
Lup 
0 0 ae 
ie : 
0 0 Gane? 
L=R,i. 


ise 





Torque Equation 
sh = Wa i, a Wa ly 


Zip, AG), 
~ @,, dt 





ace 


Tu =M(Wa)t, - MW) 


deem _ ee See m 
Ops 
T sp 
Ji 2 7 acc ny 
MB 


Structural Jacobian 


The structural jacobian for the DQO model is given by: 


= 
| 
Sys Se 
SV aoe 
Se. SS SS Ss ea SS 
eee ee ee 
om ae le ae” le o> Mag 


Jacobian Calculations 


Calculating the jacobian of the export variables with respect to the import variables 


is Straight forward with the exception of the partials with respect to the mechanical 


frequency. First of all, nothing depends on V,, hence all of its partials are zero. In the 


following derivations, the Device Jacobian is partitioned such that the voltages and 


currents are split from the mechanical speed and torque. 
I =R,CA~(B,R,V +B,5,) 


ut R,CA’'BR 
ave I vf 


- 276 - 








Calculating the partials with respect to the mechancial frequency: 


Ax =B.RV +B,5, 











R, 
= -a°(2 : -v- 24, | 











OO, “dM, 00), 
al, Ox OR, 
=R,C——+ 
dm, (tas aa oe 
where 
0 — Sp, 0.0 <0 
7A ; SP, 0 0 0 0 
00, 0 0 0 0 
0 0 0 0 0 
0 0 0 0 0 
and 
dO 
eee fC 
00),, Pp 


SB 





V Vsp 
—=M(Se)Sp, =—M(Co)Sp, 0 0 
Vip Vp 
wt "8 MC S "38 ays Sp. 0 0 
Om | Van (CoS Ps Vie Co) Pp 
0 0 0 0 
0 0 0 J 


Digi 





i MB i MB 


eee, ee ee O20 

OR, Lup lp 
ans ye cee, 7, fe) Pp 0 0 
0 0 0 0 
0 0 0 i 


The Torque equation Jacobians are given by: 
OT., Tsp OT cpu 
OV ‘Tug OV 


one = di, : OW, di, . OW, 
5 = M(Y,) st + M (i) e-MV,) 50 ML) 








OW, 

a 0 0 0 
OW, 

210 1 00 oS 
oe 
Sy 7 | 

OT irae 
ay 7h lap 


Now with respect to the mechanical frequency: 


oT, Tsp yi as =| 











JOE 00, O0,, 
ON Spee 
=—£5 
JO, O;; 
to pet. a, OY 
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Putting the Jacobian all together: 
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OWg Ox 
Bist =i 0") *{)) vad 
OW, Ox 
ane =O 7. 0:0 20) ee 
di, di. 
xa =F 10 -)0:- Bar 
a (0 J 0) ul 
00,, 00), 
Ip Ap di ee ol, 
WV_ Wo WV, I, 
Wg Ap Wp Ay 
Vy WV We 00, 
0 0 0 0 
oT, OT, oT, OT, 
~| Op Vo Wr WW, 
pipe Ore ) A) A 
WV, Wy Ve 00, 
Ap Ap Ag Alp 
WVp Vg ~ We da, 
0 0 0 0 





F-1.2 ABC Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Stator Phase A V, Gmport) I, (export) (1) Normal 

Stator Phase B V, (import) I, (export) (1) Normal 

Stator Phase C Vc (import) I. (export) (1) Normal 
Mechanical (,, (mport) T., (export) (0) Normal 

Field Voltage Vep (import) Information 

Stator Phase A Current I, (export) Information 

Stator Phase B Current Ip; (export) Information 

Stator Phase C Current I (export) Information 

Field Current I, (export) Information 


The import x;,,. and export x,,, vectors are defined by: 


ee ly 
V; Io 
Ve is 
| Ven a lr 

| ©, =e ihe 
lp 

lor 

le 


Parameters 

All Parameters are identical to the DQ0 Model 
States 

All States are identical to the DQO Model 


Equations 
Constant Definitions 


All Constants definitions are identical to the DQ0 Model 
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Angle Calculations 


For this model, the angle is the actual rotor angle of the synchronous machine: 


oe { onp,at ae 


©=S0,p, + Ov 


Note: When calculating ©,, it would be wise to limit its range to —. 


Variable Rotation and Scaling 


To convert from the V vector to the v vector, Parks transformation should be used 





M (cos(©)) 

Weg) ~Min@) 
* 3V up i 
0 


V4 
V, 
Ve 
Vep 


ofe-2) 
Hole-8) 


& 


V= 


Nl 


) 


veRV 


Solving the electrical dynamical equations 


DQO model. 
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nfofox)) 0 
“Mort )) 9 


] 
5 0 
3V 
0 wees 
2V 


The electrical dynamical equations are solved in exactly the same way as for the 





Calculating Explicit Variables 


The only difference for calculating the explicit variables are the following 


matrices: 


% 
a 
a a ee 


M (cos(©)) — M(sin(©)) 1 O 

wfeofo-28)) -u(sefo-28)) 1 
R= 20 | 21 

I<, | M| cos Shae — M| sin Se let oe) 

Isp 

0 0 0 ie 


Structural Jacobian 


The structural jacobian for the ABC model is given by: 


2 
ee ee eee ve 
z2z2z2zez2a2 
ee ee, ee 
eee ee 
Sa ee, eee ee 


Jacobian Calculations 


The only differences for calculating the jacobian matrices are the following: 


JO 


aa 
00,, Pp 
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—M(sin(O))S p, -u{ snfo-2) sp, | snf+2) bp, 





OR, 2Vsp 2n 21 
Se] nonin, -ofo-¥)fo, -ofo{o- jn 
0 0 0 
0 0 0 


—M (sin(©))S p, —M(cos(©))S p, Oe 


| oe 21 
asf oedeel lode. 6 
OW,, Isp {snl © N65, -u{ co o+22|}sp, 0 0 


0 0 0 0 





Putting the Jacobian all together: 


uu, ow, aw ah 
OV, Ws, We We dO, 
Ol, lp Ae Ap Ap 
OV, OW, We Ve Wy 
al. do Oe ao ale 
ELA AA We OV; 00, 
oT, of, of, oT, OT, 
OV, Wz We Wr dO, 
‘jm oy weecinn elem ecih 
eV, Vv, OV. OV; od, 
ol, bole Ole) foln= Col. 
OV OV, We WV, IW, 
di dle + Ole ela Oe 
WV, W, We W; 2e,| 
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F-2 Voltage Regulator Model 


This 1s a simple voltage regulator model. The voltage regulator 1s assumed to be of a 
PI type controller. This design does not have any clipping on the output waveform to 
ensure the field voltage is kept within a reasonable range. This model is intended for single 
generator operation since it has no provision for reactive power sharing with paralleled 


generators. 


F-2.1 DQ0 Model 


Interface Variables 


Terminal Potential Variable Flow Variables Type 

Line Direct Voltage Vp (amport) Information 
Line Quadrature Voltage V, (import) Information 
Reference Voltage Vf (import) Information 
Field Voltage Veep (export) Information 


The impon x,,,. and export x,,, vectors are defined by: 


Vy = [Ven 
Ximp = Vo 
Vf 
Parameters 
jie Integrating factor (1 / sec) 
ke Proportional factor (PU) 
Knoo Voltage Magnitude Conversion factor (PU) 
States 
Vyas Field Voltage 
Equations 


Calculate the terminal voltage: 


Calculate the error voltage: 


SONU 








Calculate the Field Voltage: 


Vip = | kV dt +k V.,, + Vis 


err P err 


Structural Jacobian 


The structural jacobian for the DQO model is given by: 


Jns=(N N ie 
Jacobian 

WV rp 

——=kS+k] 

Vie | 

WV ep 

Ai =e ae 


Mp _ Vey OV, 














WV, OV, W> 
Vip _ Wen OV, 
WV, OV, Wo 


Note the partials of V, with respect to Vp and Vg must be determined from the square root 


function: 
ae 
wv, PV, 
ov, Vo 
ee er 
GVige <a 


The device jacobian is given by: 


[M2 Wen Wem 
me sVn 10love 
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F-2.2 ABC Model 


Interface Variables 


Terminal Potential Variable Flow Variables Type 

Phase A Voltage V, Gmport) Information 
Phase B Voltage V, (import) Information 
Phase C Voltage V- (mport) Information 
Reference Voltage V,.¢ (import) Information 
Field Voltage Vip (export) Information 


The import X;n) and export x,,, vectors are given by: 


Ve A op = [Vin] 
ee s 
imp V. 
Ve 
Parameters 
ie Integrating factor (1 / sec) 
kp Proportional factor (PU) 
Kisc Voltage Magnitude conversion factor (PU) 
States 
V yas Field Voltage 
Equations 


Calculate and subtract out the DC offset of the common potential: 


SV al 
aes 05 
Vian = ee VG 
V50 = Va —Vo 
Veo =Ve- Vo 
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Calculate the Terminal voltage’: 


Ve = K Bc 





2 
a (Vio “Vio t Veo “Veo t+ Veo + Veo) 
Calculate the Error voltage: 


Calculate the Field Voltage: 


err 


v= fav dt +k V.,, + Vias 


Note 1: Derivation of Terminal Voltage: 


Assume phase voltages are balanced three phase: 


Vig = V-cos(O) 
2 
a= Ve cof 6 + = 


2 
Vena Ve cof © = - 


Z 2 2 Vr 
Vig Ven Veg * (1 +cos(20) + 


1 +cos(20) co “ —sin(20) in *) + 
1 + cos(20) cos = ) + sin(20) inf = } 


3 
Ve, a Ve0 te Ue. > 7 Vr 


Structural Jacobian 


The structural jacobian for the ABC model is given by: 


Jns=(N N N L} 
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Jacobian 


Vp 

= =k S +k] 
OV, pth, 
aNWis 

we —k,S ~k,I 


WVrp WV ep OV, 
OV, OV, WW, 








eo _ Vr BV, 
OV, OV, OVs 





WV _ Veo MV, 
Ve OV, We 





Note the partials of V, with respect to V,, Vs, and V- must be determined from the square 


root function: 


tt ae 
Va oe VN 
ov, _, Veo. f2 
Ven “Vv, V3 
OV coe 2 
OV FeV “NG 
De) Ce) Am ENG 


WV, 30Vig 30Ve9 30Veo 











Vv, 10, 20Y, 1 dy, 
OVs 30Vsg 30Ve9 30V eo 
OV ie ala 2 ov, 
ee SOV Slo. SOV 
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The device jacobian 1s given by: 


yee Wp WVep WVep Wep 
PANNOV, FOVn Ve aV,y 











- 289 - 






os 


+) eri ai saddoraty ot 
* 3," 
ios we - oa 
= Ante ie ial ‘ 
' bi: = rd “NE 







F-3 Prime Mover 


This is a rather crude model of a PI controller on a prime mover. The dynamics of the 


controller are assumed to dominate the response of the prime mover. 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Mechanical (,, (import) T,, (export) (0) Normal 
Information ,.¢ (import) Infonnation 


The import x,,,, and export x,,, vectors are defined by: 


exp 


a Xap =IT,,] 
Ximp = 

O,<f 

Parameters 

k, Integrating Torque factor (1 / sec) 
kp Proportional Torque factor (PU) 
W;; base frequency (rad/sec) 

Psp base System Power (watts) 

es base Machine Power (watts) 
States 

Dike mechanical torque 

Equations 


P kp 
T. =F [og ©, )dt + — (Og 0.) | 


Structural Jacobian 


The structural Jacobian for the Pinme Mover model is given by: 


fenmaen We) 
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Jacobian 


Cine ok 
JO, f PP sB ps 








ay Om OT,, 


IDn Oger 


me op Camere)! Be 
0=| San “a 











F-4 Three Phase Switch 
F-4.1 DQO Model 


Interface Variables 


Terminal Potential Variable Flow Variables 
TD1 Vp, Umport) Ty, (Export) 
TQ1 Vo, (import) Io, (Export) 
TO] V,, (Import) I,, (Export) 
TD2 Vp, (import) En. (Export) 
TQ2 Vo. (import) Io2 (Export) 
TO2 V,. (import) [42 (Export) 
SW Sy (import) 


All Interface variables are on a Per Unit (PU) Basis 


The import x;,,, and export x,,, vectors are given by: 


Vo) 


1 


] 


ty 


3 
o Oo an o 


tN 


3. 


- 291 - 


(KCL Group) Type 


(1) Normal 
(2) Normal 
(3) Normal 
(1) Normal 
(2) Normal 
(3) Normal 


Information 





Parameters 


GC... On Conductance (PU) 
Go Off Conductance (PU) 
States 


There are no states for this model. 
Equations 


The equations for the switch are very simple. First, we define the conductance G of 


the waveform: 


If Sy > 0 
Then G= Ge 
Else GG. f 


Now Generate the Conductance Matrix Gp: 


M(G) 0 0 M(-G) 0 0 0 

0 M(G) 0 0 M(-G) 0 0 

Peed ee 0 M(G) 0 0 M(-G) 0 
>| M(-G) 0 0 M(G) 0 0 0 
0 M(-G) 0 0 M(G) 0 0 

0 0 M(-G) 0 0 M(G) 0 


The export variables are simply: 


Xexp = D* imp 


Structural Jacobian 


The structural jacobian is given by: 
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N 0 0 N O O N 

ON 0 O0O N ON 

7.-}9 9 N 0 ON N 

PIN 0 0 N 0 0N 

ON O0O 0 N ON 

lo 0 N 0O0N N 

Jacobian Calculations 
The jacobian matrix is very similar to Gp: 

dG 
M(G) 90 0 M(-G) 0 0 Mp) ~Vo)5— 
a" 
dG 
0 MG) 0 0 MCG) 0 MWg ~Vo2)5—- 
a 
dG 
0 0 MG) 0 0 MCG) Mo. Vo.) 5 
- 
1 aG 
GS) 0 0 M(G) 0 0 baa Pt) ae 
: 
dG 
0 M(-G) v 0 M(G) 0 MVo2— Von) 5 
ip 
dG 
0 0 MCG) 0 0 MG) MW_-Vo)se- 
im 


Where = is determined by the If-Then-Else function. 
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F-4.2 ABC Model 


Interface Variables 


Temninal Potential Variable Flow Variables (KCL Group) Type 
TAI V,, (mport) I,, (Export) (1) Normal 

TBl V,, (Import) Iz, (Export) (2) Normal 

Tel V-, (Import) I-, (Export) (3) Normal 

TA2 V4. (Import) fd  CEXDOTL) (1) Normal 

TB2 Vz. (Import) I, (Export) (2) Normal 

Te Vc, (Import) I.) (Export) (3) Normal 

SW Sy (Import) Information 


The import x,,, and export x,,, vectors are given by: 


Va Wy, 
Vey Fy 
7 a Voy x - lo, 

Veo Ip 
Voo Lo 

Parameters 

Ge. On Conductance (PU) 

Gog Off Conductance (PU) 

States 


There are no states for this model. 


Equations 
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The equations for the switch are very simple. First, we define the conductance G of 


the waveform: 


If So) 
Then G= G,,, 
Else Gr: Gr 


Now Generate the Conductance Matnx Gp: 


M(G) 0 0 M(-G) 0 0 0 

0 M(G) 0 0 M(-G) 0 0 
Peele 0 M(G) 0 0 M(-G) 0 
>| M(-G) 0 0 M(G) 0 0 0 
0 M(-G) 0 0 M(G) 0 O 

0 0 M(-G) 0 0 M(G) 0 


The export variables are simply: 


Nexp = D* imp 
Structural Jacobian 
N 0 0 N O O WN 
0 N O0O DO N ON 
J 0 0 N QO ON N 
Nae “ON 205-0 NV 
0 N QO DO N ON 
nr 0 N 0 0NN 
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Jacobian Calculations 


The jacobian matrix is very similar to Gp: 


M(G) 0 0 M(-G) 0 0 MV, ~V.) 5 

dG 

4’ 

JG 

A aG 
JG 

dG 

0 0 M(-G) 0 0 MG) MWe-Ve))5—- 

“ 


Where — is determined by the If-Then-Else function. 
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F-5 Transmission Line 
DQ0 Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Tp, Vp, Umport) I, (Export) (1) Normal 
To Vo, (Import) Io, (Export) (2) Normal 
lie V,, (Import) I,, (Export) (3) Normal 
lee V5, Import) Tp, (Export) (1) Normal 
Tq Vo. (Import) Io2 (Export) (2) Normal 
VS V,, (Import) Ig. (Export) (3) Normal 
The import x;,,, and export x,,, vectors are given by: 
Vo) Ip, 
Vor Toy 
x _ Vor x = Io, 
‘a V2 “agle! D3 
Voo Io 
Yo Hal oo 
Parameters 
R Resistance (ohms) 
X Reactance (ohms) 
States 


(There are no states for this model) 
Equations 


Constant Definitions 


BDO 





Calculate the Export Variables 


G —Y 0 —G Y 0 
Y G 0 -Y -G 0 
0 z 0 0 at 
7 r r 
Tac yO, GeO 
-Y -G QO Y G 0 
I I 
0 0 —-— 0 0 - 
r r 
Xerp =I DXimp 
Structural Jacobian 
The Structural Jacobian is given by: 
DD OD D OO 
DD O D D OO 
Pant nO 2250 0.2 
oe kOe 002 +9720 
DD OD D OO 
0 OD 10 eK0 D| 


Jacobian Calculations 


The matrix J, is the Jacobian matrix. 
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ABC Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Wg Vo (impor) ij (Export) (1) Normal 
4 V;, (Import) I,, (Export) (2) Normal 
1s YUmporm) Ic, (Export) (3) Normal 
eS V,, (Import) Ps (Expo) (1) Normal 
Tp V,, (Import) I, (Export) (2) Normal 
a Vo, Import) coe Dei) (3) Normal 
The import x;,,, and export x,,, vectors are given by: 
ae Tay 
Ve, Ip; 
oni Vor ae ie 
imp Ve exp i 
Veo ly 
Veo | Ele 
Parameters 
R Series Resistance (ohms) 


G Parallel Conductance (mhos) 
L Inductance (hennies) 

States 

ine Phase A Inductor Current 
|e Phase B Inductor Current 

Io, Phase C Inductor Current 
Equations 


Each of the three phases can be treated independently of one another. In the 


equations which follow replace a subscripted X with the appropnate phase letter: 


First write the equations describing the phase: 


R i = 0 07) “xy 
S m1 +4 S V,.1=0 
I -=-GiVy}| |0 Gt -I|| “F7 
L jf oe 
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; oF 0 F 0 O ON Vx 
bait rab : 
eal, CS lly a te Wit ae 


Manipulating the first matrix equation, we can get an expression for the terminal 1 


ca(Seenoe}[$s0) -($-0) 


Vy 
ly, = Gc. Vy> 


ly 


current: 


Once this is known, the other variables are easy to calculate: 


Vin = Vy, RL, 


I Vy 
In = +G Wyy +3- + GV x3 
he =—ly, 


Structural Jacobian 


The Structural Jacobian is given by: 


‘. 


ae TS 2 

So. tS 
Sone fee 
Oo C-t ore 
<M a IS c= a = a = 
moonhono 
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Jacobian Calculations 


The Jacobian Matnx can be directly constructed from the first element of G,,: 


R aes 
G,, -(£s +RG 1) baa 


G Of TEC I 210 0 

(crn nO ae er eae 

0 On Ges “0 Re 
PalG 6 Oise 3G AO 0 

(UNG, <0 ef SG", 

0 Oa <IGe 10 Op Ge | 
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F-6 Constant Impedance Loads 
DQ0 Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
ce V, (amport) I, (export) (0) Normal 
To V, (amport) I, (export) (0) Normal 
dhs 9 (mport) I, (export) (0) Nornnal 


The import x;,,, and export x,,, vectors are given by: 


Vie Ls 
Mee —| V 0 bee = Io 
Vo ly 
Parameters 
R Load resistance (ohms) 
De Load reactance (ohms) 
Gond Zero Sequence conductance to ground. 
States 


(There are no states for this model) 
Equations 


First calculate the admitance 


R 

G =—— 

| ae, Ge 

ave 

R? +X? 

Now calculate the Admitance Matrix: 
G -Y () 
G,=|Y G0 
() () Ce 
X exp = Gee 
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Structural Jacobian 


The Structural Jacobian is given by: 


J ps = 


Se 
eS oo 
Rae ea 


Jacobian Calculations 


The jacobian matrix is the G,, matrix. 
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ABC Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
AW V, (amport) I, (export) (0) Normal 
16 V, Gmport) I, (export) (0) Normal 
“We V- (import) I- (export) (0) Normal 
The import x;,,, and export x,,, vectors are given by: 
WA I, 
Soe = Ve, bee —|/] B 

Ve Ic 
Parameters 
R Series Resistance (ohms) 
G Parallel Conductance (mhos) 
ig Inductance (henries) 
|e Resistance of center to Ground (ohms) 
States 
oe Phase A Inductor Current 
ee Phase B Inductor Current 
ipa Phase C Inductor Current 
Equations 


This load model can be considered to be a transmission line where all the terminals of one 
side are connected together to a resistor going to ground. As such, we can use some of 


the derivations from the transmission line model: 
R ~1 
G3, = [Es +RG 7 


S 
OF =6,{2+6] 


Ty = Gy Vy — Gy Vyo + Gish yr 


Vy 7 CU, +1, +1.)Renp 





This can be rewritten by defining the following matrix: 


7 +RenpG,) RexpGy, RenpG 
Gisc 7 RenpGir ( + RenpG;) RenpGur 
RenpG RexpGy ( + RenpG;,) 


Which allows the following equation to be written: 


if Gry O oy | G, O Oils 
Gee pa0. eG ORV AO Ga Ona, 
ie 7 0 Gal: e* 500s Geel a 
Or if we rewrite the equation: 
1p Geo Oia Ga UU Or Lana 
GeO sO a OUeV Ge oO Gi. Ober 
ie O OO Gy, Ve Oe «05 Gal ay 


The inductor current for phase X can be determined from: 


Vi = Vy — Rly 


I V 
ie =f £46 Wont B+ GReoll eyes sp 


Structural Jacobian 


The structural jacobian is given by: 


3052 





Jacobian Calculations 


The Jacobian Matrix is given by: 


Gi, 


pe Onno 


0 
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F-7 Reduction Gear 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
S1 (), (export) T, (export) (0) Normal 
$2 @, (import) T, (amport) (0) Normal 


The potential variables are measured in radians/second while the flow variables are 


measured in Newton-meters. 


The import x,,,, and export x,,, vectors are given by: 


Pe @, Be 0, 
imp ~~ es apo fie 


Parameters 

iy Number of teeth on shaft 1. 

nN, Number of teeth on shaft 2. 

nN Efficiency of Reduction Gears. 
States 


There are no states for this model. 
Equations 


The rotational speed of the shafts are proportional to the gear ratio: 


On 





The transmitted torque however, must be scaled by the efficiency. Which side of the 


equation the efficiency applies depends on the direction of the power flow: 


if T,@, > 0 
then n, 
7,=-— 
1 a 2 
else ny 
eee 
NaN 


The first zero crossing of the power should be passed back to the system as a 


suggested recalculation time. 
Structural Jacobian 


The structural jacobian is given by: 


D O 

Jos = i ” 
Jacobian Calculations 
The jacobian is given by: 
KG 
ee ny 

Poe Vol: 

dM, dT» 


The partials of 7, depend on the partial derivatives of the if-then-else function. If 


the direction of the power flow remains constant over the interval, then the partials are 


given by: 
ap 
a) 
dM, 
oti i 
OT, Tho 
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F-8 Propeller 


The relationship between the torque, angular speed, forward velocity and forces on a 
propeller are highly complex and nonlinear. While much information is known about the 
steady-state operation of propellers traversing in the forward direction, little information is 
available for nonstandard operating conditions. The classical approach is to generate K, vs J 


and , vs J curves where: 


Ee 


V 
~ nD 


F=K,J)pD n- 
Da Kp Dn 


where the variables are described in the following sections. 


The classical approach works well when V, (speed of propeller with respect to the 
fluid) and n (RPM of shaft) are both positive and n is large enough to bring J (advance 
coefficient) below about 1.5. Outside of this range, little data 1s provided for most 
propellers. The classical approach breaks down completely when the shaft speed is zero 
and J is infinite. Furthermore, there is no way to differentiate between backing down (n 
and V, both negative) and having forward way on (n and V, both positive). The method 
used for this model is better suited for simulation studies because it essentially uses the 
angle of attack on the propeller blade as the argument for the thrust and torque coefficients. 
This model is based on work conducted at the Naval Ship Research and Development 
Center, Annapolis, MD by D. W. Baker and C. L. Patterson and reflects data and theory 
developed by I. Ya. Miniovich. 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Mechanical (,, (import) T,,, (€Xport) (0) Normal 
Hydrodynamic u (import) F (export) (0) Normal 


Note: units are radians/second, Newton-meters, meters/second, and Newtons 


The Import X,,,, and Export X,,, vectors are defined by: 
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Parameters 

D Diameter of Prop (meters) 

Ww Wake Fraction (PU) 

fe) Density of water (kg/m’) 

Ci) Thrust Coefficient matrix (unlimited rows by 2 columns) 
first column is © in radians [-n 1] 
second column is Thrust Coefficient in PU. 

CoQ) Torque Coefficient matrix (unlimitd rows by 2 columns) 
first column is © in radians [-n 1] 
second column is Torque Coefficient in PU. 

States 


(There are no states for this model) 


Equations 


@ = atan2(nD,V,) 
F =-C,(@)pD°(V, +n°D’) 
T,, = Co(©)pD*(V, +n°D*) 


Normally, C;{) and C)() are specified as data points. Hence some type of 
interpolations scheme is needed to determine the value of theses functions at intermediate 


points, as well as the value of the first derivative. 
Device Structural Jacobian 


The Device Structural Jacobian for all waveforms is given by: 
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N WN 
a 
. i 4 
Jacobian Calculations 


Calculate partials with respect to intermediate variables 


OF 








as (1 = yy 
Tm 
Of ep 
Oi lade 
oe non 
Ty _1 
dw, 2m On 
Calculate the partials 
OF dC;(©) 
BES =-pD (2c; (O)V, +(V> +n°D*)———— av, | 
OF 2 2 2 0C;(O) 
oF = pb 2Cx(0)D n+(Vi+n'D’) >, | 





oe) 


OT, 
— 2Co(O)V, +(V, +n°D*) 


ifs aC6) 
== pD {2 (Q)D*n + (V2 +n°D?)—=— | 


Calculate the partial of © with respect to V, and n. 


sOlee a 
ne 
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V 
P 


00 nD 


on | +(% f 
nD 
Of course, the partials of C(O) and Co(O) with respect to © must be determined 
from the interpolation scheme used. 


Putting all of this together: 








oT,, OT, 
ne OW, ou 
- | OF OF 
OW, ou 
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F-9 Ship Dynamics Model 


Interface Variables 

Terminal Potential Variable Flow Vanables (KCL Group) Type 

Ship Hydrodynamics u (import) F (export) (0) Normal 
Velocity u is measured in meters/second while force F is measured in Newtons. 


The import x;,,, and export x,,, vectors are given by: 


imp 

X imp oz [u ] Xexp i LF] 
Parameters 
fe) Density of Salt Water (kg/m’) 1025.9 kg/m’ @ 15° C. 
v Kinematic Viscosity of Water (m*/sec) 1.19x10° m*/sec @15° C. 
G Acceleration of Gravity (m/sec’) 9.80665 m/sec’ 
L Length of Ship (m) 
As Surface Area of Ship (m’) 
m Mass of Ship (kg) 
M wad Added Mass Multiplier (PU) (normally between 1.0 and 1.10) 
Cc. Correlation Allowance 
CAR.) Matrix of Frictional Drag Coefficients (2Xn) or (3xn) 

column 1 are Reynolds Number Values 

column 2 are Frictional Drag Coefficient Values 

column 3 are optional first derivative values of the curve 

Note: Values should be provided for negative Reynolds Numbers 
CAF,) Matrix of Residual Drag Coefficients (2x) or (3xn) 

column 1 are Froude Number Values 

column 2 are Residual Drag Coefficient Values 

column 3 are optional first derivative values of the curve 

Note: Values should be provided for negative Froude Numbers 
States 


There are no States associated with this device 


- 313 - 





Tae 


q 


7 | 








Equations 


The basic equations are given by: 


fF = 


we fere 


Cr= CAR) +C(F,)+C, 





: d 
ja =P wAgCy +m, - 


The only potential difficulty is performing the evaluation of the drag coefficients 
CAR,) and C,(F,). 


Structural Jacobian 
The structural jacobian is given by: 
Jos = IN] 


Jacobian Calculations 


0 L_ ( dCAR.) 1 dC,(F,) 
Jp = 5A) M(u)M (u)| ~M| — ‘Te AF + 2M (u)M(C(R,)) | + 


-1 
ins 








dC4R,) dC, (F,) 


where S is the integration matrix and —,— (—,;—) must be determined by either 


differentiating the CAR.) (C,(F,)) curve or by interpolating the third colum if so provided. 


ae 





F-10 Pulse Generator 
Interface Variables 


Terminal Potential Variable Flow Variables Type 
VO V (export) Information 


The import x;,,, and export x,,, vectors are given by: 


X imp = {] Xexp a [V] 
Parameters 
Vug Value of V when off 
V on Value of V when on 
tp Matrix of pulse times (2x77,): 
column | are pulse on times. 
column 2 are pulse off times. 
unlimited number of n, rows. 
States 
None 
Equations 


If the time f falls between an on and an off time then V = V,,,, otherwise V = Voy. 


If a discontinuity falls within the time interval, the earliest discontinuity time should 
be passed back as a recommended recalculation time. The time of the next discontiuity 


after the time interval should also be made available to the system solver. 
Structural Jacobian 

There is no structural jacobian matrix for this device. 
Jacobian Calculations 


There 1s no jacobian matrix for this device. 





F-11 Induction Motor 
F-11.1 ABC Model 


Interface Variables 


Terminal Potential Variable Flow Variables (KCL Group) Type 
Phase A V,, (Import) I, (Export) (1) Normal 
Phase B V, (Import) I, (Export) (1) Normal 
Phase C V- Umport) I, (Export) (1) Normal 
Neutral V, (Import) I, (Export) (1) Normal 
Mechanical w,, (Import) T,, (Export) (0) Normal 


Voltages are in volts, currents in amps, angle speed in radians per second and 


Torque in Newton-meters. 


The import x; and export x,,, vectors are given by: 


Va I, 
V, lp 
Ximp = Vo 5 Io 
Vo I 
L®, Tm, 

Parameters 

R; Stator Resistance (ohms) 

Rp Rotor Resistance (reflected to Stator) (ohms) 

rs Stator Reactance (ohms) 

AF Rotor Reactance (reflected to Stator) (ohms) 

GP Mutual Reactance (ohms) 

J Moment of Inertia (Kg-m’) 

‘wo Base Frequency (radians per second) 

Pp Pole Pairs 

B Windage Torque Factor (Newton-Meters-second) 
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States 
S) Electrical Rotor Angle (radians) 
Equations 


First calculate the electrical angle: 
i= 0,+ | p,@adt 


Now specify the following stator voltage and current vectors: 


WE. I, 
Vs =| Va—Vo ls =| 4p 
Ve V, Ic 


The Rotor voltages and currents as reflected on the stator are: 


4 , 

Var Lap 

, , - , 
Ve =| Ver ig =| /or 
, 4 

Vee lor 


Calculate the inductances: 





xX 

7 b=! 
Ops Ws 
X 

Vey L. ==M 
(,. 3 


The Stator induction matrix is given by: 


1 ] 
peek =e ==) 
( ls a) 2 ms 9; ms 
1 ] 
L, = 5b ms (Let) — 5 Ems 
] ] 
a=] Soil web Lege 
9) ms D ms ( ls aid 


The Rotor induction matrix is given by: 
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] l 
eee IE fey, 
( lr ey, a ms 2 ms 


; I ; ] 
Lp a = 5 ems (L, a is _ 9 je 
L ae (L,°+L_.) 
9 ms a; ms lr ms | 


The Mutual induction matrix is given by: 


M (cos(8))L,,, M cos tr < ))e M (cos 0 — 2 a 


ee, Sie (cof = = )} M (cos(®))L,,, M cos 0 + = a 


ihe (co < = }}bn M (co = 2 \z M (cos(®))L,,, 


The Rotor and Stator resistance matrices are given by: 


I 0 0 10 0 
R’=R,0 1 0 R,=R{0 I 0 
CaO ey 007 


r 


The system of equations which need to be solved can be expressed as: 
ale Rese Solow @ 7S lees a ail 
Ve = Sc ((heae R’ +S"'L,’ i 
For a squirrel cage induction motor vg. = 0. Hence the rotor curents can be solved 


for: 


ip = Aggls 
, —]} a | \T 
Ba a a Sie) 
Arps = —(S ee a5 ae (aoe 


Now the stator currents can be solved for: 
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Aaah ee Sale Ate 
Vom Ave 
The neutral current is solved using KCL; 
i=, 15-1, 


The electrical torque is given by: 


2X yP Tgp Ter’ 
T, =-—— | M(,)M] 1,,’-—--— |+ 
: 3, { Us) 1, 2e- w2 


Loe 
MU) Iy’—"— J 
I I 
CLM (ly 6] sn + 
3 ? ? / ? ? , 
‘/ tM LM or —Iep )+MUp)M cr Lip )AMU)M Lin — Lp )) cos 
The full torque equation is given by: 


T, =-T.+(JS"+Blqo, 


Structural Jacobian 


z2z222a 
z=z22a 
Se 
ee ee 
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Jacobian Calculations 
The jacobian elements corresponding to the electrical variables are easy to 
calculate. If we partition the jacobian matrix as follows: 
Ass AssU, J ew 
Jp = Ur Ass Ur AssU, OF ier 
Ite J reU, Sow 


where 
—] 
U,=| -1 
—] 


The remaining matrices Uw, Jz, and J7,) are not as easy to calculate: 





Med AGE 
Take the partial wrt @,,: 
(A dis 
je = ( ss) | a 
Jo, ° S 90,, 
Oi O(Lp’A 
ee eee Bese “a As Aas) 
JO,, O,, 


Lisp Ars a —Lep (SR, a joe (aoe 








O(Lsp’A O(Lsp’Y AL sp’ 
aos) = A Lg (SR, +L’) i 5p, (Re thn’ Ee’) 


M (sin(9))L,,,5 Pp o{sin(¢ P Jens, oso = 2 }}nsp, 


arr 21 , ; 27 
——_ = Msn sa 3 } ens, M(sin(9))L,,.5 P, a sin(¢ aE 3. ms Pp 


Msn +2 } ens, o sn(o — 2 Jens, M (sin(8))L,,,.5 P, 
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Appendix G: WAVESIM Program Files 


WAVESIM consists of ten program source code (.c) files plus three header (.h) files 
and must be linked to the standard math hbrary. Six of the program source code files and 
two of the header files are specific to WAWVESIM while the remaining files contain 


application independent code. 


WAVESIM SPECIFIC FILES 
wavesim.c Main Executive Routine 
wavesim.h Definition of WAVESIM Structures 
waveinit.c Initialization of Structures 


Read device.def 


waveinit.h Define Initial values of all system parameters 
waveread.c Read and interpret Input File 

wavebld.c Build System 

wavewrit.c Write MATLAB .M file 

wavewrta.c Write MATLAB .M file (continued) 


APPLICATION INDEPENDENT FILES 


ioliba.c String Manipulation Routines 
iolbia.h declarations of ioliba.c routines 
get _file.c Prompt for and open files 
getline.c Obtain string from an input stream 
filebase.c File name manipulation routines 
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G-1 Main Program File: wavesim.c 


wavesim.c contains the main routine which performs the executive functions for 


WAVESIM: 
1. Initialize device definitions 
2. Print the Header 
3. Open Files 
4. Read Input File 
5. Build the System 
6. Write The Output File 
7. Close the Files 
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init devices 


read file 
build system 


write file 





G-2 System Initialization: waveinit.c 


waveinit.c contains the following routines for initializing the system: 


init devices 


read device def 

print _system_base 

print device def 

print matrix 

print structural _ jacobian 


read terminal 


read parameter 


read st ate 


read f unction 


read structural jacobian 


strip white 


read matrix 
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Sets Default Values 

Calls read _ device def to read in 
device.def file. 

Debug handler. 

Reads in device. def file 

Prints system base parameters 

Prints a Device Definition 

Prints a Matrix 


Prints a Structural Jacobian 


Interprets TERMINAL — subordinate 
command 


Interprets PARAMETER subordinate 
command 


Interprets STATE subordinate 


command 


Interprets FUNCTION subordinate 
command 


Interprets STRUCTURAL 
JACOBIAN subordinate command 


Strips all blanks, tabs, retums from a 
string. 


Reads a matrix. 





The hierarchy for the routines in waveinit .c is given by: 


init devices 
read device def 
read terminal 
read parameter 
read matrix 
read state 
read function 
read structural jacobian 
strip white 
print system base 
print_device_ def 
Present ematrix 


print structural jacobian 
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G-3 Reading Input File: waveread.c 


waveread.c contains the following routines for reading in an input file: 


read file 


read f ile device 


read f 1 le default 


read file node 


EeaG.s lest ime 


read file debug 


print _ debug 
print time 
print device 
print system 
print _node 


finish reading file 


2325= 


Controls other routines for reading 
input files 

Debug handler. 

Reads and _  =Interprets device 
command from input file. 


Reads and Interprets default 
command from input file. 


Reads and Interprets node command 
from input file. 


Reads and Interprets time command 
from input file. 


Reads and __ Interprets 
command from input file. 


debug 


Print debug flag status. 

Print simulation time parameters. 
Print device characteristics. 

Print system characteristics. 

Print node characteristics. 

Generate cross references and 


generally finish the process. of 
developing the initial system. 





The hierarchy for the routines in waveread.c 1s given by: 


read file 

read file device 

print device 
print _matrix 

read matrix 

read file default 
print _system_base 

read file node 
print node 

read file time 
print time 

read file debug 
print debug 

fonireh reading f2le 

print_system 
print _device 


print _node 
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G-4 Building the System: wavebld.c 


wavebld.c contains the following routines for building the system: 


build system 


build system identify 


build system xref 


build system structural _ jacobian 


build system blocks 


find block 


print system identify 


sj_add 


sj sub 


print system block_sJjac 


print block 
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Identifies System variables and 
equations 

Builds cross references 

Builds System Structural Jacobian 


Reduces system 


Identifies system variables and 


equations. 


Builds cross references within the 


system. 


Builds the system _ structural 


jacobian matrix. 


Identifies the sequence of blocks for 


solving system. 


Attempts to find a block of a given 


Size. 
Prints system information. 


Adds two Structural Jacobian 


codes. 


Subtracts two Structural Jacobian 


codes. 


Prints the block owner of each 
element in the system sructural 


jacobian matrix. 


Prints information about a block. 





The hierarchy for the routines In wavebld.c is given by: 


build system 
build system_identify 
print system identify 
build system xref 
build system_structural jacobian 
sj _add 
sj sub 
print structural jacobian 
build system blocks 
find block 
print block 


print structural jacobian 
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G-5 Writing MATLAB M-File: wavewrit.c and wavewrta.c 


wavewrit.c and wavewrta.c contain the following routines for writing the 


output MATLAB M-File: 


write file 


write file header 


WELceura le sinictaalize 


write file time loop 


~ 


write file plot_variables 


write fi remLOOLer 


wr ite f1 le so lve block 


Calls routines to 


MATLAB M-File. 


other generate 


Prints header information to 


MATLAB M-File. 


Prints System Initialization 


parameters to MATLAB M-File. 
Pnnts Time Loop algorithm to 
MATLAB M-Fue. 


Prints algorithm 
variables to MATLAB M-File 


to plot system 


Footer information to 


MATLAB M-File. 


Prints 


Prints algorithm for solving block to 
MATLAB M-File. 


The hierarchy for the routines in wavewrit.c and wavewrta.c is given by: 


write file 
write f Tlevheader 
Wetter nile wntagt ize 


write file time loop 


write file solve block 


write file footer 
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G-6 Application Independent Files 


Several Support files containing special C functions are required by WAVESIM. 
These files were written by the author independently of WAVESIM and may contain 
routines unused by WAVESIM. 


G-6.] ioliba.c 


i0oliba.c contains a number of functions for manipulating strings. The functions 
used by WAVESIM are: 


Stoda Converts a string to an array of double precision floating 
numbers. 

strextract Extracts the nth word of a string 

strsplit Returns the remainder of a string after the nth word. 

strstrip Strips a string of leading and trailing spaces, tabs, and 


Carriage returns 


strncmpa Case insensitive version of strncemp for comparing the 


first n characters of two Strings. 


strcmpa Case insensitive version of strcmp for comparing two 


strings. 


- 330 - 





G-6.2 getline.c 


getline.c contains functions for reading in lines from a file and automatically 


implementing the following features: 


g 
Ms, 


>. 


Comment Lines beginning with #, !, or % are ignored. 
Blank Lines are ignored. 
Lines can be continued with ... or \. 


If the first word of the line is INCLUDE followed by a filename, lines are read 
from that file. (NOTE: A check for recursive INCLUDE statements is 
performed to prevent infinite loops) 


Carriage Returns are truncated. 


The functions used by WAVESIM are: 


init _strm Initialization function 


get_line Function used to read the lines in. 
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G-6.3 get_file.c 


get _file.c contains functions for opening input and output files. The following 


features are implemented: 


l. 


Two strings are passed: a default filename string and an argument filename 


string. Either or both strings may be empty. 


If the argument filename is specified, the functions attempt to open that file. 
If opening that file is unsuccessful, the user is prompted to enter a new 


filename. 


If the argument filename is empty, the user is prompted to enter a filename. If 
the default filename is not empty, it 1s offered to the user as the default name 
of the file. If opening the file entered by the user is unsuccessful, the user is 


prompted to enter a new filename. 


Users can exit the routine without opening a file by entering only a q when 


prompted for a filename 


Users can obtain a directory listing by entering a ? followed by whatever file 


specification the user desires. 
Leading and trailing spaces in filenames are truncated. 


Whatever filename 1s successfully opened is passed back as a new default 


filename. 


The function used by WAVESIM 1s: 


get input file Function for opening an Input File 


G-6.4 filebase.c 


filebase.c contains functions for stripping an extension off of a filename and for 


returning the extension of a filename. The following features are implemented: 


ie 


An extension is defined as all the characters after the last period (.) found after 
the last directory delimiter (\ for IBM-DOS and / for UNIX). If no such period 


is found, the extension does not exist. 


The function used by WAVESIM is: 


extract base Extract the base filename (without extension) 
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G-7 Makefile 


The UNIX make utility greatly eases the task of developing programs by only compiling 


those files which have changed since the last compilation. Here is the Makefile used to 
generate WAVESIM on a VAXstation 3100: 


# Makefile for wavesim 

# 

# For Revision 2.0 

# 

FILES = wavesim.c ioliba. 
waveinit.c waveread. 
wavebld.c 

OBJ = wavesSim.o ioliba. 
waveinit.o waveread. 


c getline.c get file.c filebase.c \ 
Cc wavewrit.c wavewrta.c \ 


© -getiane 0 get file.o filebases,o ~\ 
O wavewrit.o wavewrta.o \ 


wavebld.o 
HEADER= wavesim.h ioliba.h 


CFLAG = -g 
SOMPILE = Ce 
# 
# 
# 
wavesim: $(HEADER) $ (OBJ) 

S(COMPILE) -o wavesim $(CFLAG) $(OBJ) -lm 
wavesim.o: $(HEADER) wavesim.c 

S (COMPILE) -c $(CFLAG) wavesim.c 
TOLioa.-O. si1Oliba.~h- tOliba.c 

S (COMPILE) -c $(CFLAG) ioliba.c 
getline.o: getline.c 

S$(COMPILE) -c $(CFLAG) getline.c 


Beemer le .o. ge. filewc 

S$ (COMPILE) -c $(CFLAG) 
filebase.o: filebase.c 

S (COMPILE) -c $(CFLAG) 
waveinit.o: $ (HEADER) 

S$ (COMPILE) -c $(CFLAG) waveinit.c 
waveread.o: $(HEADER) waveread.c 

S$ (COMPILE) -c $(CFLAG) waveread.c 
wavewrit.o: $S(HEADER) wavewrit.c 

S$ (COMPILE) -c $(CFLAG) wavewrit.c 
wavewrta.o: S$(HEADER) wavewrta.c 

S (COMPILE) -c $(CFLAG) wavewrta.c 
wavebld.o: $(HEADER) wavebld.c 

$ (COMPILE) -c $(CFLAG) wavebld.c 


Gev vtec. ¢ 


filebase.c 


# 
i 
# 
al 9 yom 

jane So (P LLES) 


oo. 


waveinit.h waveinit.c 
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